Georges Aazan
email: georges.aazan@lsv.fr

Antoine Girard
email: antoine.girard@l2s.centralesupelec.fr

Luca Greco
email: luca.greco@l2s.centralesupelec.fr

Paolo Mason
email: paolo.mason@l2s.centralesupelec.fr

Stability of Shuffled Switched Linear Systems: A Joint Spectral Radius Approach

Keywords: Switched systems, Shuffled stability, Joint spectral radius, Lyapunov methods, Automata theoretic techniques

A switching signal for a switched system is said to be shuffled if all modes of the system are activated infinitely often. In this paper, we develop tools to analyze stability properties of discrete-time switched linear systems driven by shuffled switching signals. We introduce the new notion of shuffled joint spectral radius (SJSR), which intuitively measures how much the state of the system contracts each time the signal shuffles (i.e. each time all modes have been activated). We show how this notion relates to stability properties of the associated switched systems. In particular, we show that some switched systems that are unstable for arbitrary switching signals can be stabilized by using switching signals that shuffle sufficiently fast and that the SJSR allows us to derive an expression of the minimal shuffling rate required to stabilize the system. We then present several approaches to compute lower and upper bounds of the SJSR using tools such as the classical joint spectral radius, Lyapunov functions and finite state automata. Several tightness results of the bounds are established. Finally, numerical experiments are presented to illustrate the main results of the paper.

Introduction

Switched systems are dynamical systems with several modes of operation, each mode being described by a differential (continuous-time) or difference (discrete-time) equation. At each instant, the active mode is determined by a so-called switching signal. Switched systems are very useful in applications to describe in a faithful manner the execution of control algorithms on distributed computing infrastructures and thus for taking into account the constraints related to the sharing of computing and communication resources [START_REF] Alur | Compositional modeling and analysis of multi-hop control networks[END_REF][START_REF] Donkers | Stability analysis of networked control systems using a switched linear systems approach[END_REF][START_REF] Su | Stability of a class of linear switching systems with applications to two consensus problems[END_REF].

The question of the stability of switched systems has been extensively studied for almost three decades. Early works have mostly focused on stability for switching signals that are arbitrary or that satisfy some minimum or average dwell-time condition (see e.g. [START_REF] Liberzon | Basic problems in stability and design of switched systems[END_REF][START_REF] Liberzon | Switching in systems and control[END_REF][START_REF] Sun | Switched linear systems: control and design[END_REF][START_REF] Lin | Stability and stabilizability of switched linear systems: a survey of recent results[END_REF] and the references therein). For the class of discretetime switched linear systems with arbitrary switching signals, a powerful notion for analyzing stability is that of joint spectral radius (JSR) [START_REF] Jungers | The joint spectral radius: theory and applications[END_REF], which allows to characterize the worst-case asymptotic growth of infinite products of matrices. However, computing the JSR is, in general, a difficult task [START_REF] Tsitsiklis | The Lyapunov exponent and joint spectral radius of pairs of matrices are hardwhen not impossible -to compute and to approximate[END_REF], and several methods have been proposed to compute tight lower and upper bounds [START_REF] Gripenberg | Computing the joint spectral radius[END_REF][START_REF] Blondel | Computationally efficient approximations of the joint spectral radius[END_REF][START_REF] Vankeerberghen | JSR: A toolbox to compute the joint spectral radius[END_REF]. In the past decade, a significant amount of research has been carried out to analyze the stability of switched systems subject to constrained switching. Constraints on the switching signal are usually described by finite-state automata and methods for handling such constraints through Lyapunov functions have been presented in several works [START_REF] Lee | Uniformly stabilizing sets of switching sequences for switched linear systems[END_REF][START_REF] Kozyakin | The Berger-Wang formula for the Markovian joint spectral radius[END_REF][START_REF] Philippe | Stability of discrete-time switching systems with constrained switching sequences[END_REF]. An alternative approach to analyze stability with such constrained switching signals is by computing the so-called constrained joint spectral radius [START_REF] Dai | A Gel'fand-type spectral radius formula and stability of linear constrained switching systems[END_REF][START_REF] Kozyakin | The Berger-Wang formula for the Markovian joint spectral radius[END_REF][START_REF] Xu | Approximation of the constrained joint spectral radius via algebraic lifting[END_REF].

There are however classes of constrained switching signals that cannot be described using classical finite-state Preprint submitted to Automatica 2 May 2022

automata. An example of such a class is that of shuffled switching signals. A switching signal is said to be shuffled if each mode of the switched system is activated an infinite number of times. It is known from formal language theory [START_REF] Baier | Principles of model checking[END_REF] that a non-trivial set of sequences cannot satisfy both liveness and safety properties. Then, the set of shuffled switching signals, which satisfies the liveness property, cannot be characterized using classical finite state automata that can only describe safety properties. To the best of our knowledge, shuffled stability has first been studied in [START_REF] Kozyakin | Algebraic unsolvability of problem of absolute stability of desynchronized systems[END_REF] where it is proved that the set of stable shuffled switched linear systems is not semi-algebraic. However, no further characterization of shuffled stability is provided in that paper. A necessary and sufficient condition can be found in [START_REF] Gurvits | Stability of discrete linear inclusion[END_REF] in terms of the maximal spectral radius of an infinite set of matrices. However, its computation is non-trivial and not discussed in that paper, which therefore does not provide any practical way to check shuffled stability. In [START_REF] Wang | Stability analysis of switched linear systems defined by regular languages[END_REF], it is shown that a notion of robust shuffled stability is equivalent to stability for arbitrary switching signals. In [START_REF] Girard | Lyapunov functions for shuffle asymptotic stability of discrete-time switched systems[END_REF], a Lyapunov characterization of shuffled stability and a converse Lyapunov result have been established. While the latter paper provides a method to check shuffled stability of switched linear systems, it does not provide tight estimates of the convergence rate in particular with respect to the characteristics of the shuffling phenomena. Providing such estimates is the main contribution of this paper.

Stability analysis of shuffled switched systems is also interesting from the point of view of applications. An example of such applications is multi-agent consensus with switching communication topologies. Consider a symmetric and connected communication graph where only one edge is active at each time instant and is selected by a switching signal. If one uses a shuffled switching signal, then each edge is activated infinitely often so that the union of future communication graphs is connected at all time, and it is well-known (see e.g. [START_REF] Moreau | Stability of multiagent systems with timedependent communication links[END_REF][START_REF] Blondel | Convergence in multiagent coordination, consensus, and flocking[END_REF]) that the consensus is asymptotically reached. On the contrary, if some edges are never activated the union of future communication graphs may be disconnected and the consensus is not attainable. Another potential application is the design of observers for switched systems where the dynamics in each mode is unobservable. Clearly, the dynamics of the estimation error cannot be stable for arbitrary switching signals, since keeping the same mode activated all the time makes the system unobservable (see e.g. [START_REF] Sun | Switched linear systems: control and design[END_REF]). However, considering shuffled switching signals, it may be possible to design asymptotically convergent observers [START_REF] Tanwani | Observability for switched linear systems: characterization and observer design[END_REF]. This type of observer design problems is natural when considering a system with distributed sensors communication over a shared network. More generally, shuffled switching signals are of interest to describe applications where multiple components rely on a shared resource and where the access to this resource must be granted to each of the component infinitely often. An example of such application is shown in an illustrative case study at the end of this paper. An additional motivation for considering shuffled switching signals is that they constitute a simple example of more general ω-regular languages [START_REF] Baier | Principles of model checking[END_REF] that can be described using either non-deterministic Büchi automata, deterministic Rabin automata, deterministic Muller automata (which are finite state automata with various acceptance conditions), or using Linear Temporal Logic formulas. Moreover, it has been shown in [START_REF] Wang | Stability analysis of switched linear systems defined by regular languages[END_REF] that the stability of a switched linear system whose switching signal belongs to an arbitrary ω-regular language is actually equivalent to the shuffled stability of a lifted switched linear system. Being able to analyze the stability of such systems gives the possibility to deal with complex applications involving control and scheduling co-design [START_REF] Weiss | Automata based interfaces for control and scheduling[END_REF][START_REF] Alur | Regular specifications of resource requirements for embedded control software[END_REF].

In this paper, we provide a new tool to analyze shuffled stability of discrete-time switched linear systems. The main contributions of this paper are as follows. We introduce a notion of shuffled joint spectral radius (SJSR), which intuitively measures how much the state of the system contracts each time the signal shuffles (i.e. each time all modes have been activated). We establish several properties of the SJSR and show how it relates to stability properties of switched systems driven by shuffled switching signals. In particular, we show that some switched systems that are unstable for arbitrary switching signals can be stabilized by using switching signals that shuffle sufficiently fast, the minimal shuffling rate being related to both the JSR and the SJSR. We also present two approaches to compute approximation of the SJSR. The first approach is based on the JSR of a finite set of matrices and allows us to compute asymptotically tight lower and upper bounds. The second approach is based on Lyapunov functions and automata theoretic techniques and allows us to compute upper bounds. The Lyapunov approach draws inspiration from [START_REF] Girard | Lyapunov functions for shuffle asymptotic stability of discrete-time switched systems[END_REF] but we consider a different underlying automata which allow us to prove the tightness of our upper bound. Our Lyapunov conditions also resemble those of [START_REF] Kundu | On stability of discretetime switched systems[END_REF], however, in that work, the relation of such Lyapunov functions to shuffled switching signals was not established.

The organization of the paper is as follows. Section 2 presents the necessary background on switched systems and shuffled switching signals. In Section 3, the SJSR is introduced and several of its properties are established. Section 4 shows the relationship between the SJSR and stability properties of a switched system driven by shuffled switching signals. In Section 5, we present two approaches to compute approximations of the SJSR. Finally, two numerical examples are used in Section 6 to illustrate the main results of the paper.

Notations: R + 0 denotes the set of non-negative real numbers. I n ∈ R n×n denotes the identity matrix. We use • to denote an arbitrary norm on R n and the associated induced matrix norm defined for M ∈ R n×n by M = sup

x =0 M x
x . The spectral radius of a matrix

M is denoted ρ(M ), note that ρ(M ) = lim k→+∞ M k 1/k . Given a sequence of matrices (M k ) k∈N , with M k ∈ R n×n , we define for K 1 , K 2 ∈ N, with K 1 ≤ K 2 K2 k=K1 M k = M K2 • • • M K1 , if K 2 > K 1 , M K1 , if K 2 = K 1 .
Given a set of matrices M ⊆ R n×n and α ∈ R, we define αM = {αM | M ∈ M}. Given a set I, we use 2 I to denote its powerset, that is the set of subsets of I, including the empty set ∅ and I itself.

Preliminaries

Let us consider a discrete-time switched system of the following form:

x(t + 1) = A θ(t) x(t), t ∈ N (1)
where x : N → R n is the trajectory and θ : N → I is the switching signal belonging to the class S(I) of arbitrary switching signals. I = {1, . . . , m}, with m ≥ 2, is the finite set of modes and A = {A i ∈ R n×n |i ∈ I} is a collection of matrices indexed by the modes. For a switching signal θ, let A θ,0 = I n , and

A θ,T = T -1 t=0 A θ(t) , ∀T ≥ 1.
Given an initial state x 0 ∈ R n , the trajectory defined by (1) with x(0) = x 0 is unique and is denoted x(., x 0 , θ), it satisfies for all t ∈ N, x(t, x 0 , θ) = A θ,t x 0 .

In this paper, we focus on a particular class of constrained switching signals called shuffled, i.e. signals for which each mode in I is activated infinitely often. A formal definition is given in [START_REF] Girard | Lyapunov functions for shuffle asymptotic stability of discrete-time switched systems[END_REF] as: Definition 1. A switching signal θ :

N → I is shuffled if ∀i ∈ I, ∀T ∈ N, ∃t ≥ T : θ(t) = i.
The sequence of shuffling instants (τ θ k ) k∈N is defined by τ θ 0 = 0 and for all k ∈ N,

τ θ k+1 = min t > τ θ k ∀i ∈ I, ∃s ∈ N : τ θ k ≤ s < t and θ(s) = i .
The shuffling index κ θ : N → N is given by

κ θ (t) = max{k ∈ N| τ θ k ≤ t}. t θ(t) t κ θ (t) τ θ 0 τ θ 1 τ θ 2 τ θ 3
Fig. 1. A switching signal (with m = 2) and the associated shuffling instants and shuffling index [START_REF] Girard | Lyapunov functions for shuffle asymptotic stability of discrete-time switched systems[END_REF].

The shuffling rate γ θ is defined as

γ θ = lim inf t→+∞ κ θ (t) t .
The set of shuffled switching signals is denoted S s (I).

Intuitively, the period between two successive shuffling instants τ θ k and τ θ k+1 is a period of minimal duration over which all modes are activated at least once. The shuffling index κ θ (t) counts the number of times the signal has shuffled (i.e. the number of shuffling instants) up to time t. The shuffling rate γ θ characterizes the frequency of shuffling instants over all time. Note that for all k ∈ N, τ θ k ≥ km and therefore, for all t ∈ N, we have 0 ≤ κ θ (t) ≤ t m , from which it follows that γ θ ∈ [0, 1 m ]. An example of shuffled switching signal is shown in Figure 1.

It is well known (see e.g. [16, Corollary 1.1]) that the system (1) is asymptotically stable for arbitrary switching signals if and only if ρ(A) < 1, where ρ(A) is the joint spectral radius (JSR) of A, formally defined for a bounded set of matrices M ⊆ R n×n by:

ρ(M) = lim k→+∞   sup    k j=1 M j 1/k M j ∈ M, j = 1, . . . , k      .
However, the stability condition based on the JSR turns out to be conservative in the case of constrained switching sequences. Some works have adapted the notion of joint spectral radius to switching signals generated by deterministic [START_REF] Dai | A Gel'fand-type spectral radius formula and stability of linear constrained switching systems[END_REF][START_REF] Kozyakin | The Berger-Wang formula for the Markovian joint spectral radius[END_REF][START_REF] Xu | Approximation of the constrained joint spectral radius via algebraic lifting[END_REF] or stochastic [START_REF] Gurvits | Stability of discrete linear inclusion[END_REF][START_REF] Chitour | On the gap between deterministic and probabilistic joint spectral radii for discrete-time linear systems[END_REF] finite state automata. However, none of these constraints allow to deal with shuffled switching signals. Therefore, we introduce, in the next section, a new notion of joint spectral radius adapted for shuffled switching signals.

Shuffled Joint Spectral Radius

In this section, we introduce the shuffled joint spectral radius and establish several of its properties. Definition 2. Let ρ > ρ(A), the Shuffled Joint Spectral Radius relative to (A, ρ) (ρ-SJSR for short) is defined as

λ(A, ρ) = lim sup k→+∞   sup θ∈Ss(I) A θ,τ θ k ρ τ θ k 1/k   . (2)
To better understand the characteristics of the ρ-SJSR, we are going to analyze some properties of the function λ(A, •). We first show that the ρ-SJSR of A can be seen as the JSR of an infinite but bounded set of matrices. For ρ > ρ(A), let us consider the following set of matrices:

M ρ =        1 ρ N N k=1 A j k , j 1 , . . . , j N ∈ I, N ∈ N, ∀i ∈ I, ∃k, j k = i,
and ∀k = N, j k = j N        .
Since N is unbounded, it is clear that M ρ generally has an infinite number of elements. The following Lemma shows that this infinite set is bounded for ρ > ρ(A). Lemma 1. For all ρ > ρ(A), M ρ is a bounded set of matrices.

Proof. Let ρ > ρ(A), then ρ(

1 ρ A) = ρ(A) ρ < 1.
Therefore, from [6, Theorem 1 (b)], we get that 1 ρ A is a left convergent product (following the terminology in [START_REF] Berger | Bounded semigroups of matrices[END_REF]). Since the elements of M ρ are products of matrices belonging to the set 1 ρ A, it follows from [6, Theorem 1 (a)] that M ρ is a bounded set of matrices.

Since M ρ is a bounded set of matrices, its JSR is well defined. Moreover, in any matrix of M ρ , all modes in I appear at least once with the last mode j N appearing only once. Hence, the set M ρ coincides with the set of all possible matrices

1 ρ τ θ 1 A θ,τ θ 1 for θ ∈ S s (I). Lemma 2. For all ρ > ρ(A), λ(A, ρ) = ρ(M ρ ). Proof. Let ρ > ρ(A) and k ∈ N with k ≥ 1, observe that for every θ ∈ S s (I) the sequence (θ j ) j∈N of switching signals in S s (I) defined by θ j (t) = θ(t + τ θ j ) for t, j ∈ N is such that A θ,τ θ k = k-1 j=0 A θj ,τ θ j 1 and τ θ k = k-1 j=0 τ θj 1 . (3)
On the other hand, given a sequence (θ j ) j∈N in S s (I) one can construct a switching signal θ ∈ S s (I) such that

τ θ i = i-1 j=0 τ θj 1 for i ≥ 1 and θ(t) = θ i (t -τ θ i ) for τ θ i ≤ t < τ θ i+1
, so that (3) holds true. Then, it follows that sup θ∈Ss(I)

A θ,τ θ k ρ τ θ k 1/k = sup θ0,...,θ k-1 ∈Ss(I)    k-1 j=0 A θj ,τ θ j 1 ρ τ θ 0 1 +•••+τ θ k-1 1    1/k = sup M1,...,M k ∈Mρ k j=1 M j 1/k .
Then, taking the limit superior on both sides when k goes to infinity yields the expected result.

Due to Lemma 2, the ρ-SJSR inherits several properties of the JSR.

Proposition 1. For all ρ > ρ(A), λ(A, ρ) = lim k→+∞   sup θ∈Ss(I) A θ,τ θ k ρ τ θ k 1/k   . (4)
Moreover, the value λ(A, ρ) does not change if we replace the induced matrix norm in [START_REF] Alur | Compositional modeling and analysis of multi-hop control networks[END_REF] or in (4) by any other matrix norm.

Proof. The result is a direct consequence of Lemma 2 and of the properties of the JSR stated in [16, Lemma 1.2] and [START_REF] Jungers | The joint spectral radius: theory and applications[END_REF]Page 10].

We provide further useful properties of λ(A, •) in the proposition below. Proposition 2. One of the following properties holds:

(i) The function ρ → λ(A, ρ) is decreasing, takes values in (0, 1) and, for all ρ(A) < ρ 1 ≤ ρ 2 , λ(A, ρ 2 ) ≤ ρ 1 ρ 2 m λ(A, ρ 1 ). ( 5
)
(ii) For all ρ > ρ(A), λ(A, ρ) = 0 and for all θ ∈ S s (I) A θ,t = 0, for all t ≥ τ θ n .

Proof. Let ρ > ρ(A), similar to the proof of Lemma 1, it follows from [6, Theorem 1] applied to the set of matrices

1 ρ A that there exists C 1 ≥ 1 such that 0 ≤ A θ,T ρ T ≤ C 1 , ∀θ ∈ S(I), ∀T ∈ N. (6)
By considering T = τ θ k , raising to the power 1/k and taking the supremum over shuffled switching signals, we find 0 ≤ sup

θ∈Ss(I) A θ,τ θ k ρ τ θ k 1/k ≤ C 1/k 1 .
Taking the limit of all terms yields λ(A, ρ) ∈ [0, 1].

Let us assume that there exists ρ > ρ(A) such that λ(A, ρ ) = 0 and let us consider an arbitrary shuffled switching signal θ ∈ S s (I). Let us consider the sequence (θ j ) j∈N of switching signals in S s (I), defined as in the proof of Lemma 2. From Lemma 2, ρ(M ρ ) = 0 and therefore the joint spectral radius of any finite subfamily of M ρ is also zero. In particular, we get for the following subfamily of n elements:

ρ 1 ρ τ θ 0 1 A θ,τ θ 0 1 , . . . , 1
ρ τ θ n-1 1 A θ,τ θ n-1 1 = 0 .
By applying [16, Proposition 2.1] and using (3) we get A θ,τ θ n = 0. It follows that A θ,t = 0, for all t ≥ τ θ n , for all θ ∈ S s (I). Then, we get that λ(A, ρ) = 0 for all ρ > ρ(A).

Now let us assume that for all

ρ > ρ(A), λ(A, ρ) > 0. Let ρ(A) < ρ 1 ≤ ρ 2 . By recalling that τ θ k ≥ km,
we have for all k ∈ N and for all θ ∈ S s (I)

A θ,τ θ k ρ τ θ k 2 = A θ,τ θ k ρ τ θ k 1 ρ 1 ρ 2 τ θ k ≤ A θ,τ θ k ρ τ θ k 1 ρ 1 ρ 2 km .
Raising to the power 1/k and taking the supremum over shuffled switching signals yields

sup θ∈Ss(I) A θ,τ θ k ρ τ θ k 2 1/k ≤ sup θ∈Ss(I) A θ,τ θ k ρ τ θ k 1 1/k ρ 1 ρ 2 m .
Now we take the limit of both terms and we get [START_REF] Barabanov | Lyapunov indicator of discrete inclusions. 1. 2. 3[END_REF], which implies that ρ → λ(A, ρ) is decreasing.

Then, let us assume that there exists ρ 2 > ρ(A) such that λ(A, ρ 2 ) = 1. It follows from ( 5) that for all

ρ 1 ∈ (ρ(A), ρ 2 ), λ(A, ρ 2 ) < λ(A, ρ 1 ), which contradicts the fact that λ(A, ρ 1 ) ∈ [0, 1]. Hence, λ(A, ρ) ∈ (0, 1) for all ρ > ρ(A).
In order to get rid of the dependence of ρ in the ρ-SJSR, it is natural to introduce the following definition. Definition 3. The Shuffled Joint Spectral Radius (SJSR) of A is defined as

λ(A) = lim ρ→ρ(A) + λ(A, ρ). (7)
Since by Proposition 2, λ(A, •) is bounded and nonincreasing in ρ, the right limit at ρ(A) in ( 7) exists and the SJSR is well-defined. We can now show some properties of the SJSR. Proposition 3. The SJSR enjoys the following properties:

(i) λ(A) belongs to [0, 1] and is independent of the choice of the norm; (ii) For all K ∈ R, K = 0, we have λ(KA) = λ(A).

Proof. The first statement follows from [START_REF] Blondel | Convergence in multiagent coordination, consensus, and flocking[END_REF] and by the properties of λ(A, •) proved in Propositions 1 and 2.

Concerning the second item, let K ∈ R, K = 0, we have by ( 7)

λ(KA) = lim ρ→ρ(KA) + λ(KA, ρ) = lim ρ→ρ(A) + λ(KA, |K|ρ)
where the second equality comes from the property of the JSR, ρ(KA) = |K|ρ(A), see e.g. [16, Proposition 1.2]. Furthermore, from (2), one can deduce that for all

ρ > ρ(A), λ(KA, |K|ρ) = λ(A, ρ), therefore λ(KA) = λ(A).
Let us remark that while the JSR ρ(A) belongs to R + 0 , the SJSR λ(A) is always a value in [0, 1]. Intuitively, while the JSR provides an estimate of the contraction rate (when ρ(A) < 1) or of the expansion rate (when ρ(A) > 1) of the system state at each time step for arbitrary switching signals, the SJSR measures how much additional contraction of the state is obtained each time the signal shuffles. Theorem 1 in the next section provides theoretical ground to this interpretation.

Shuffled switched systems and ρ-SJSR

In this section, we show how the SJSR relates to stability properties of switched linear systems with shuffled switching signals. From Proposition 2, we already know that if λ(A, ρ) = 0 for some ρ > ρ(A) then all trajectories of (1) will stay at 0 after n shuffling instants. So, we focus in this section on the case when λ(A, ρ) > 0 for all ρ > ρ(A). In particular, we show how the ρ-SJSR allows us to compute bounds on the shuffling rate ensuring stability. The next result clarifies the relationship between the ρ-SJSR and the behavior of the trajectories of (1). Theorem 1. For all ρ > ρ(A), for all λ ∈ (λ(A, ρ), 1], there exists C ≥ 1 such that

x(t, x 0 , θ) ≤Cρ t λ κ θ (t) x 0 , ∀θ ∈ S s (I), ∀x 0 ∈ R n , ∀t ∈ N. (8)
Conversely, if there exists C ≥ 1, ρ ≥ 0 and λ ∈ [0, 1] such that (8) holds, then either ρ > ρ(A) and λ ≥ λ(A, ρ), or ρ = ρ(A) and λ ≥ λ(A).

Proof. We start by proving the direct result. Let ρ > ρ(A) and λ ∈ (λ(A, ρ), 1]. By definition of λ(A, ρ), there exists k 0 ≥ 1 such that sup

θ∈Ss(I) A θ,τ θ k ρ τ θ k 1/k ≤ λ, ∀k ≥ k 0 .
It follows that

A θ,τ θ k ≤ ρ τ θ k λ k , ∀θ ∈ S s (I), ∀k ≥ k 0 . (9)
Then, let C 1 ≥ 1 be such that (6) holds. In particular, for T = τ θ k and for shuffling switching signals, we obtain from (6) that

A θ,τ θ k ≤ C 1 ρ τ θ k , ∀θ ∈ S s (I), ∀k ∈ N. (10)
Then, let C 2 = C 1 λ -k0 , it follows from ( 9) and ( 10) that

A θ,τ θ k ≤ C 2 ρ τ θ k λ k , ∀θ ∈ S s (I), ∀k ∈ N. (11)
Let θ ∈ S s (I), t ∈ N, and k = κ θ (t), we have

A θ,t = A θ ,t-τ θ k A θ,τ θ k
where θ ∈ S s (I) is given by θ (s) = θ(τ θ k + s), for all s ∈ N. By (6), we get that

A θ ,t-τ θ k ≤ C 1 ρ t-τ θ k . (12)
Then, let C = C 1 C 2 , by [START_REF] Dai | A Gel'fand-type spectral radius formula and stability of linear constrained switching systems[END_REF] and ( 12), we get

A θ,t ≤ A θ ,t-τ θ k A θ,τ θ k ≤ Cλ κ θ (t) ρ t .
Hence, (8) holds.

We now prove the converse result. By definition of induced matrix norm, ( 8) is equivalent to the following:

A θ,t ≤ Cρ t λ κ θ (t) , ∀θ ∈ S s (I), ∀t ∈ N. (13)
Since λ ∈ [0, 1] and κ θ (t) ∈ N we also have

A θ,t ≤ Cρ t , ∀θ ∈ S s (I), ∀t ∈ N.
Raising the previous terms to the power 1/t and taking the supremum over shuffling switching signals yields sup θ∈Ss(I)

A θ,t 1/t ≤ C 1/t ρ, ∀t ∈ N.
Using the fact that for every θ ∈ S(I) and all t ∈ N there exists θ t ∈ S s (I) that coincides with θ up to time t, the equation above can be re-written as sup θ∈S(I)

A θ,t 1/t ≤ C 1/t ρ, ∀t ∈ N.
Taking the limit as t goes to infinity, one obtains that ρ(A) ≤ ρ. Recalling (13) and fixing t = τ θ k , we have

A θ,τ θ k ≤ Cρ τ θ k λ k , ∀θ ∈ S s (I), ∀k ∈ N which is equivalent to A θ,τ θ k ρ τ θ k 1/k ≤ C 1/k λ, ∀θ ∈ S s (I), ∀k ∈ N. ( 14
)
If ρ > ρ(A), taking the supremum over all shuffled switching signals and the limit as k goes to infinity yields

λ(A, ρ) ≤ λ. If ρ = ρ(A), then for all ρ > ρ(A), (14) gives
A θ,τ θ k ρ τ θ k 1/k ≤ C 1/k λ, ∀θ ∈ S s (I), ∀k ∈ N.
Then, it follows that for all ρ > ρ(A), λ(A, ρ ) ≤ λ and therefore λ(A) ≤ λ.

The previous theorem provides a bound on the growth of the state and can be used to derive conditions for stabilization of system (1) using shuffled switching signals with a minimal shuffling rate. Corollary 1. Assume λ(A, ρ) > 0 for every ρ > ρ(A).

Let θ ∈ S s (I), if there exists ρ > ρ(A) such that γ θ > -ln(ρ) ln(λ(A,ρ)) , then

lim t→+∞ x(t, x 0 , θ) = 0, ∀x 0 ∈ R n . ( 15
)
Proof. If λ(A, ρ) > 0 and γ θ > -ln(ρ) ln(λ(A,ρ)) , then there exists λ ∈ (λ(A, ρ), 1) such that γ θ > -ln(ρ) ln(λ) and > 0 such that < γ θ + ln(ρ) ln(λ) . Then, from the definition of shuffling rate, there exists t 0 ∈ N such that

κ θ (t) t > γ θ -2 , ∀t ≥ t 0 .
From Theorem 1, there exists C ≥ 1 such that for all x 0 ∈ R n , for all t ≥ t 0 .

x(t, x 0 , θ)

≤ Cρ t λ κ θ (t) x 0 ≤ Cρ t λ (γ θ -2 )t x 0 ≤ Cρ t λ (-ln(ρ) ln(λ) + 2 )t x 0 = Cλ 2 t x 0 ,
from which (15) follows.

A practical consequence of the previous corollary is that even if the switched system is unstable for arbitrary switching (i.e. if ρ(A) > 1), it can be stabilized by shuffling sufficiently fast. Note that since for all θ ∈ S s (I), the shuffling rate γ θ ≤ 1 m , it follows that stabilization by shuffling can only be done if there exists ρ > ρ(A) such thatln(λ(A,ρ)) ln(ρ) > m.

Remark 1. As a consequence of Corollary 1, we have that if λ(A) ∈ (0, 1), γ θ > -ln(ρ(A)) ln(λ(A)) implies [START_REF] Gurvits | Stability of discrete linear inclusion[END_REF]. However, we are unable to show that, in general, the function ρ → -ln(ρ) ln(λ(A,ρ)) is non-decreasing, hence replacing the assumption of Corollary 1 by γ θ > -ln(ρ(A)) ln(λ(A)) may increase the conservativeness of the result.

Approximation of the ρ-SJSR

Similar to the JSR, the exact computation of the ρ-SJSR appears in general out of reach, and an interesting problem is to look for theoretical and numerical methods allowing to estimate the ρ-SJSR. Such estimates would have an impact on the applicability of the stability results, Theorem 1 and Corollary 1. From Lemma 2, it follows that λ(A, ρ) could be approximated by computing bounds on ρ(M ρ ). However, off-the-shelf algorithms for computing approximations of the JSR only apply to finite set of matrices [START_REF] Vankeerberghen | JSR: A toolbox to compute the joint spectral radius[END_REF]. A lower bound of λ(A, ρ) can still be computed by computing a lower bound of the JSR of a finite subfamily of M ρ . However, a similar approach does not allow to compute upper bounds. In the next subsection, we provide lower and upper bounds based on the computation of the joint spectral radius of a finite set of matrices. Methods to compute upper bounds for the ρ-SJSR based on Lyapunov techniques are then developed in Section 5.2.

Asymptotic estimate

In this section, we provide asymptotically tight bounds of λ(A, ρ) for a large enough ρ. We first recall the following classical result, which roughly speaking ensures, for any bounded set of matrices M, the existence of a matrix norm such that ρ(M) is approximated by the maximum norm of the matrices in M. Proposition 4 ([6, Lemma 2]). Let M be a bounded set of matrices. The following equality holds

ρ(M) = inf • max M ∈M M ,
where the infimum is taken among all matrix norms induced from norms in R n .

When it exists, an induced norm • that satisfies ρ(M) = max M ∈M M is called extremal for M. We refer to [START_REF] Jungers | The joint spectral radius: theory and applications[END_REF]Theorem 2.2], and also to [START_REF] Barabanov | Lyapunov indicator of discrete inclusions. 1. 2. 3[END_REF] and [START_REF] Kozyakin | Algebraic unsolvability of problem of absolute stability of desynchronized systems[END_REF], for sufficient conditions and for necessary and sufficient conditions for the existence of extremal norms.

Let us consider the set N I of products of matrices where all modes in I appear exactly once:

N I = m k=1 A j k j 1 , . . . , j m ∈ I,
∀i ∈ I, ∃k ∈ {1, . . . , m}, j k = i .

The following proposition provides a lower bound for the ρ-SJSR. Proposition 5. For all ρ > ρ(A)

λ(A, ρ) ≥ ρ(N I ) ρ m .
Proof. Let us remark that we have 1 ρ m N I ⊆ M ρ . Then, from the definition of the JSR, we get

ρ(M ρ ) ≥ ρ 1 ρ m N I = ρ(N I ) ρ m .
Lemma 2 then allows us to conclude.

It is natural to ask whether the lower bound provided in Proposition 5 can prove asymptotically tight. The following theorem provides an upper bound for the ρ-SJSR and answers to this question. Theorem 2. The following results hold true.

(i) For all K > ρ(N I ), there exists R ≥ ρ(A) such that for all ρ ≥ R,

λ(A, ρ) ≤ K ρ m . (16)
(ii) We have the asymptotic estimate

lim ρ→+∞ ρ m λ(A, ρ) = ρ(N I ). ( 17
)
(iii) If there exists a norm • * that is extremal for N I , then there exists R ≥ ρ(A) such that for all ρ ≥ R,

λ(A, ρ) = ρ(N I ) ρ m . ( 18
)
Proof. Let us fix K > ρ(N I ). By Proposition 4 there exists an induced matrix norm • * such that M * ≤ K for every M ∈ N I . Considering now the set Mρ = ρ m M ρ , we have N I ⊆ Mρ . We want to show that M * ≤ K for every M ∈ Mρ \ N I , provided that ρ is large enough. For this purpose, we first take any R > ρ(A). Then, again by Proposition 4, there exists an induced matrix norm • * * such that A * * ≤ R for every A ∈ A. Taking K ≥ 1 such that A * ≤ K A * * for every matrix A ∈ R n×n we obtain

A θ,i * ≤ K A θ,i * * ≤ K Ri , ∀θ ∈ S(I), ∀i ∈ N.
Observe now that every element of Mρ \ N I takes the form

1 ρ τ θ 1 -m A θ,τ θ 1 for some θ ∈ S s (I) such that τ θ 1 > m. Setting R = max R, K K Rm+1 we have for every ρ ≥ R 1 ρ τ θ 1 -m A θ,τ θ 1 * ≤ 1 ρ τ θ 1 -m K Rτ θ 1 = K Rm R ρ τ θ 1 -m ≤ K Rm R R τ θ 1 -m ≤ K Rm+1 R ≤ K.
This proves that M * ≤ K for every M ∈ Mρ \ N I , hence for every M ∈ Mρ , if ρ ≥ R. We deduce that

λ(A, ρ) = 1 ρ m ρ( Mρ ) ≤ 1 ρ m sup M ∈ Mρ M * ≤ K ρ m ,
for ρ ≥ R, which proves [START_REF] Jungers | The joint spectral radius: theory and applications[END_REF].

Setting K = ρ(N I ) + for an arbitrary > 0 we then get the existence of R ≥ ρ(A) such that, for all ρ ≥ R,

0 ≤ ρ m λ(A, ρ) -ρ(N I ) ≤ ,
where the inequality on the left follows from Propositions 5. Letting tend to zero, this implies [START_REF] Kozyakin | Algebraic unsolvability of problem of absolute stability of desynchronized systems[END_REF].

Finally, by taking • * extremal for N I in the argument above, we obtain that (16) holds true with K = ρ(N I ). Together with Proposition 5, this implies [START_REF] Kozyakin | The Berger-Wang formula for the Markovian joint spectral radius[END_REF].

In the case where all the matrices A i are invertible, a practical criterion to estimate the lower bound R in Theorem 2 is given by the following lemma. Before stating the lemma, we define, for each i ∈ I, the set N i of products of N matrices in A with N < m, where each mode in I appears at most once and i appears exactly once:

N i = N k=1
A j k j 1 , . . . , j N ∈ I, N < m, ∀k = k , j k = j k , and ∃k, j k = i .

Lemma 3. Assume all the matrices in A are invertible. Let K ≥ ρ(N I ) and let • * be an induced matrix norm such that M * ≤ K for every M ∈ N I . Then, setting

R = max i∈I,M ∈Ni M A i M -1 * , (19)
we have that R ≥ ρ(A) and (16) (and (18), if K = ρ(N I )) holds true for ρ ≥ R.

Proof. Following the proof of Theorem 2, it is enough to show that A θ,τ θ 1 * ≤ KR τ θ 1 -m for every θ ∈ S s (I) such that τ θ 1 > m. The product of matrices

A θ,τ θ 1 = A θ(τ θ 1 -1) • • • A θ(0)
is shuffling once, therefore all the modes appear at least once (with θ(τ θ 1 -1) appearing exactly once). The total number of repetitions in the product is given by n = τ θ 1 -m > 0. Suppose that the first repeated mode is θ(j 1 ). In order to remove A θ(j1) , using the fact that all matrices A i are invertible, we multiply A θ,τ θ 1 on the right by

Q 1 = (A θ(j1-1) • • • A θ(0) ) -1 A -1 θ(j1) (A θ(j1-1) • • • A θ(0) ). Then, A θ,τ θ 1 Q 1 = A θ (1) ,τ θ (1)
1

where the switching signal θ (1) is given by θ (1) (t) = θ(t) until the instant j 1 -1, and with θ (1) (t) = θ(t + 1) afterwards. It is clear that

τ θ (1) 1 = τ θ 1 -1.
By repeating this process n times, we get a switching signal θ (n) without repetitions before its first shuffling instant (i.e. τ θ (n)

1 = τ θ 1 -n = m).
Then, we finally have

A θ,τ θ 1 n k=1 Q n+1-k = A θ (n) ,τ θ (n) 1 , with A θ (n) ,τ θ (n) 1 ∈ N I . It follows that A θ,τ θ 1 = A θ (n) ,τ θ (n) 1 n k=1 Q -1 k ,
where for all k = 1, . . . , n,

Q -1 k ∈ i∈I {M A i M -1 |M ∈ N i },
and therefore

Q -1 k * ≤ R, where R is given by (19). Moreover, since A i ∈ N i , we have that R ≥ max i∈I A i * ,
which implies by Proposition 4 that R ≥ ρ(A). Finally, from submultiplicativity of the induced matrix norm, we get that for all θ ∈ S s (I)

A θ,τ θ 1 * ≤ A θ (n) ,τ θ (n) 1 * n k=1 Q -1 k * ≤ KR n ,
which yields the expected result since n = τ θ 1 -m.

We have seen with Theorem 2 and Lemma 3 that, under certain conditions, we can find tight bounds on the ρ-SJSR. In particular, if one can find a norm • * that is extremal for N I , ( 19) provides an R ≥ ρ(A) such that for all ρ ≥ R the exact value of ρ-SJSR is given by [START_REF] Kozyakin | The Berger-Wang formula for the Markovian joint spectral radius[END_REF]. For all ρ ∈ [ρ(A), R) we can still use the lower bound provided by Proposition 5, but we need an effective way to compute an upper bound. To address this problem, in the next section we present two Lyapunov-based methods.

Before that, we end this section by exhibiting a class of switched systems for which the value of the ρ-SJSR can be explicitly determined for all values of ρ > ρ(A). Let us consider a set of commuting matrices A = {A 1 , • • • , A m }, it is well known (see e.g. [START_REF] Blondel | On the accuracy of the ellipsoid norm approximation of the joint spectral radius[END_REF]) that in that case the JSR is equal to the maximal value of the spectral radii of the matrices of A, that is ρ(A) = max(ρ(A 1 ), . . . ρ(A m )). Interestingly, in that case, we are also able to determine the expression of the ρ-SJSR: Proposition 6. For a finite set of commuting matrices

A = {A 1 , • • • , A m }, we have for all ρ > ρ(A), λ(A, ρ) = ρ(A 1 • • • A m ) ρ m .
Proof. Since the matrices commute, all matrices N I are equal to

A 1 • • • A m . Then, Proposition 5 provides a lower bound on λ(A, ρ) with ρ(N I ) = ρ(A 1 • • • A m ).
We need to show that it is also an upper bound.

From Proposition 4 we know that, for all ρ > ρ(A), there exists a (submultiplicative) matrix norm • * such that ρ > max i∈I ( A i * ). In the case of commuting matrices, the product A θ,τ θ k can be rearranged as follows

A θ,τ θ k = (A 1 • • • A m ) k A i1 1 • • • A im m where mk + i 1 + • • • + i m = τ θ k . Whereby we have A θ,τ θ k * ρ τ θ k ≤ A 1 • • • A m k * A 1 i1 * • • • A m im * ρ km ρ i1 • • • ρ im .
Moreover, recalling that

Aj * ρ < 1 for j = 1, . . . , m, we get A θ,τ θ k * ρ τ θ k ≤ A 1 • • • A m k * ρ km .
By raising the previous expression to the power 1 k , by taking the supremum over θ ∈ S s (I) and the limit for k going to infinity, we get

λ(A, ρ) ≤ ρ(A 1 • • • A m ) ρ m .

Computing bounds with Lyapunov techniques

This section details two methods for computing upper bounds on the JSR and the ρ-SJSR. Both approaches are based on Lyapunov functions and automata theoretic techniques. While the first approach is computationally more tractable, the second approach is shown to be tight.

Lyapunov function indexed by I

Our first method is based on a Lyapunov function indexed by the set of modes I. It can be seen as an extension of Theorem 1 in [START_REF] Girard | Lyapunov functions for shuffle asymptotic stability of discrete-time switched systems[END_REF] in order to compute bounds on the JSR and the ρ-SJSR (in [START_REF] Girard | Lyapunov functions for shuffle asymptotic stability of discrete-time switched systems[END_REF], only a stability certificate is provided).

Theorem 3. If there exist V :

I ×R n → R + 0 , α 1 , α 2 , ρ > 0 and λ ∈ [0, 1] such that the following inequalities hold true for every x ∈ R n α 1 x 2 ≤ V (i, x) ≤ α 2 x 2 , i ∈ I (20) V (i, A i x) ≤ ρ 2 V (i, x), i, i ∈ I, i = i (21) V (i + 1, A i x) ≤ ρ 2 V (i, x), i ∈ I \ {m} (22) V (1, A m x) ≤ ρ 2 λ 2m V (m, x), (23)
then the bound (8) holds.

Proof. The proof follows the main lines of that of Theorem 1 in [START_REF] Girard | Lyapunov functions for shuffle asymptotic stability of discrete-time switched systems[END_REF]. Let us consider an initial condition x 0 ∈ R n and a shuffled switching signal θ, let us denote x(•) = x(•, x 0 , θ). Let η : N → I be defined by η(0) = 1 and the following rules:

   if θ(t) = η(t), then η(t + 1) = η(t); if θ(t) = η(t) and η(t) = m, then η(t + 1) = η(t) + 1; if θ(t) = η(t) and η(t) = m, then η(t + 1) = 1.
(24) An equivalent description of the evolution of η can be given in terms of a finite state automaton, like the one shown in Figure 2 for m = 3. Inspecting the dynamics of η, one can easily see that η takes transitions from m to 1 an infinite number of times if and only if θ is a shuffled switching signal. Moreover, there is at least one such transition every m shuffling instants. Let us consider W : N → R + 0 , defined by

W (t) = V (η(t),x(t)) ρ 2t
, for all t ∈ N. It follows from ( 21), ( 22) and ( 23) that W (t + 1) ≤ W (t), thus W (t) is non-increasing. Proceeding as in the proof of Theorem 1 in [START_REF] Girard | Lyapunov functions for shuffle asymptotic stability of discrete-time switched systems[END_REF], we can show that

W (τ θ jm ) ≤ λ 2jm W (0), ∀j ∈ N. Let t ∈ [τ θ jm , τ θ (j+1)m
), then κ θ (t) ≤ (j + 1)m. From the monotonicity of W (t) we have

W (t) ≤ W (τ θ jm ) ≤ λ 2κ θ (t)-2m W (0).
From [START_REF] Lee | Uniformly stabilizing sets of switching sequences for switched linear systems[END_REF] we have W (0) ≤ α 2 x 0 2 and we get for all

t ∈ N V (η(t), x(t)) ≤ α 2 ρ 2t λ 2κ θ (t)-2m x 0 2 .
From ( 20) we also have, for all t ∈ N,

x(t) 2 ≤ V (η(t),x(t)) α1
. Whence, by taking C = λ -m α2 α1 , we obtain the bound [START_REF] Blondel | Computationally efficient approximations of the joint spectral radius[END_REF].

The previous result allows us to compute upper-bounds on the JSR and the ρ-SJSR. However, there is some conservativeness that is mostly due to the fact that the Lyapunov function is guaranteed to contract only every m shuffling instants. In the following subsection, we show an approach that makes it possible to compute a Lyapunov function that contracts at each shuffling instant. This also allows us to compute tight upper bounds on the JSR and the ρ-SJSR.

Lyapunov function indexed by 2 I

Our second approach is also based on the use of a Lyapunov function and of an automaton. The main difference is that we use an automaton that visits a given state at each shuffling instant (instead of at most each m shuffling instants as in the previous section). In order to do that, we need to use a Lyapunov function indexed on the powerset of I. While this significantly increases the complexity of the approach, 2 I having exponentially more elements than I, this allows us to propose a method that provides tight bounds on the JSR and the ρ-SJSR. Theorem 4. If there exist V : (2

I \ {I}) × R n → R + 0 , α 1 , α 2 , ρ > 0 and λ ∈ [0, 1] such that the following in- equalities hold true for every x ∈ R n α 1 x 2 ≤ V (J, x) ≤ α 2 x 2 , ∀J I (25) V (J ∪ {i}, A i x) ≤ ρ 2 V (J, x), if J ∪ {i} = I (26) V (∅, A i x) ≤ ρ 2 λ 2 V (J, x), if J ∪ {i} = I (27)
then the bound (8) holds. Conversely, if the matrices A i are invertible, for all i ∈ I and the bound (8) holds for some ρ > 0, λ ∈ [0, 1] and C ≥ 1, then there exists a function V : (2 I \ {I}) × R n → R + 0 such that the inequalities (25), ( 26) and (27) are satisfied.

Proof. We first prove the direct part. Let us consider an initial condition x 0 ∈ R n and a shuffled switching signal θ, let us denote x(•) = x(•, x 0 , θ). Let η : N → 2 I \ {I} be defined by η(0) = ∅ and the following rules:

if η(t) ∪ {θ(t)} = I, then η(t + 1) = η(t) ∪ {θ(t)}; if η(t) ∪ {θ(t)} = I, then η(t + 1) = ∅.
(28) An equivalent description of the evolution of η can be given in terms of a finite state automaton, like the one shown in Figure 3 for m = 3. It is straightforward to see from the dynamics of η that η(t) = ∅ if and only if t = τ θ k for some k ∈ N. Hence, the automaton state ∅ is visited at each shuffling instant.

Now, let us consider

W : N → R + 0 , defined by W (t) = V (η(t),x(t)) ρ 2t
, for all t ∈ N. It follows from ( 26) and ( 27) that for all t ∈ N, W (t + 1) ≤ W (t), thus W (t) is nonincreasing. Also since for all k ∈ N, we have η(

τ θ k ) = ∅, (27) gives W (τ θ k ) ≤ λ 2 W (τ θ k -1) for every k ∈ N, k ≥ 1.

The latter relation and the monotonicity of

W (t) implies that W (τ θ k ) ≤ λ 2 W (τ θ k-1 )
. Hence, by induction we conclude that

W (τ θ k ) ≤ λ 2k W (0), ∀k ∈ N. Let t ∈ [τ θ k , τ θ (k+1)
), then κ θ (t) = k, from the monotonicity of W (t) we have:

W (t) ≤ W (τ θ k ) ≤ λ 2κ θ (t) W (0).
From this point on the proof is similar to the final part of that of Theorem 3. Therefore we finally obtain the bound [START_REF] Blondel | Computationally efficient approximations of the joint spectral radius[END_REF]. We prove now the converse result. We first consider the case ρ = 1. For J ⊆ I with J = ∅, let M J consist of all finite products of N matrices in A, with N ∈ N, where each mode in J appears at least once, and the last mode j N belongs to J and appears exactly once:

∅ {1} {2} {3} {1, 2} {2, 3} {1, 3} 2 1 2, 3 1, 3 1, 2 2 1 1
M J =        N k=1 A j k , j 1 , . . . , j N ∈ I, N ∈ N, ∀i ∈ J, ∃k, j k = i, j N ∈ J and ∀k = N, j k = j N        .
In particular, M I is the set of all possible matrices A θ,τ

1 θ for θ ∈ S s (I). Since ρ = 1, it follows from (8) that for all k ∈ N A θ,τ k θ ≤ Cλ k , ∀θ ∈ S s (I) ,
which is equivalent to

k j=1 M j ≤ Cλ k , ∀M 1 , . . . , M k ∈ M I .
From [6, Lemma 2] applied to the family of products of matrices in {1 λ M | M ∈ M I }, it follows that there exists a norm • * in R n such that the corresponding induced matrix norm, also denoted by • * , satisfies sup

M ∈M I M * ≤ λ.
Then, for all subsets J I and x ∈ R n we define

V (J, x) = sup M ∈M I\J M x 2 * .
Let us first prove [START_REF] Philippe | Stability of discrete-time switching systems with constrained switching sequences[END_REF]. For J I, let M J be an arbitrary element of M I\J , then V (J, x) ≥ M J x 2 * . All matrices M J are invertible, • * is equivalent to • , and 2 I \ {I} is a finite set, therefore there exists α 1 > 0 such that M J x 2 * ≥ α 1 x 2 , for all J I, x ∈ R n . Also, it follows from (8) with ρ = 1 and λ ∈ [0, 1] that for all J I, for all M ∈ M I\J , M ≤ C. Then, • * being equivalent to • , it follows that there exists α 2 > 0 such that V (J, x) ≤ α 2 x 2 , for all J I, x ∈ R n .

To prove [START_REF] Su | Stability of a class of linear switching systems with applications to two consensus problems[END_REF], we first notice that every product M A i with M ∈ M I\(J∪{i}) belongs to I\J . Hence,

V (J ∪ {i}, A i x) = sup M ∈M I\(J∪{i}) M A i x 2 * ≤ sup M ∈M I\J M x 2 * = V (J, x).
In order to prove [START_REF] Sun | Switched linear systems: control and design[END_REF], we use the fact that if J I is such that

J ∪ {i} = I then A i ∈ M I\J = M {i} . Consequently, V (∅, A i x) = sup M ∈M I M A i x 2 * ≤ λ 2 A i x 2 * ≤ λ 2 sup M ∈M I\J M x 2 * = V (J, x).
This concludes the proof in the case ρ = 1.

Let us consider now the general case ρ = 1. For the set of matrices 1 ρ A, the corresponding dynamics satisfies (8) with ρ = 1 and we can consider the function V defined above. It directly follows from the properties shown above and the fact that V is homogeneous of degree two that ( 25)-( 26)-( 27) hold true. The proof is complete.

As a direct consequence of the previous results and Theorem 1 we have the following corollary. Corollary 2. Let ρ > 0 and λ ∈ [0, 1] such that there exists a function V satisfying the conditions of Theorem 3 or Theorem 4, then either ρ > ρ(A) and λ ≥ λ(A, ρ), or ρ = ρ(A) and λ ≥ λ(A). Conversely, if the matrices A i are invertible, for all i ∈ I, then for all ρ > ρ(A), for all λ ∈ (λ(A, ρ), 1], there exists a function V satisfying the conditions of Theorem 4.

Corollary 2 shows that upper-bounds of the JSR and of the ρ-SJSR can be found by computing Lyapunov functions satisfying the conditions in Theorems 3 and 4. Limiting the search to quadratic Lyapunov function, the conditions ( 20)-( 21)-( 22)-( 23) or ( 25)-( 26)- [START_REF] Sun | Switched linear systems: control and design[END_REF] can straightforwardly be translated to linear matrix inequalities (LMIs) for which efficient solvers exist. However, the tightness of the conditions in Theorem 4 is lost when constraining the Lyapunov functions to be quadratic.

Numerical examples

In this section, we illustrate the main results of the paper with two numerical examples 1 .

The first example shows an application of Proposition 5 to compute a lower bound of the ρ-SJSR, and of Theorem 2, Theorems 3 and 4 to compute upper bounds of the ρ-SJSR. The second example illustrates the result of Corollary 1 where a switched system that is unstable for arbitrary switching is stabilized by shuffling sufficiently fast. The switched systems considered in both examples are adapted from [START_REF] Girard | Lyapunov functions for shuffle asymptotic stability of discrete-time switched systems[END_REF].

Computing bounds of the ρ-SJSR

Let us consider a switched system in R 3 with 2 modes, where A = {A 1 , A 2 } with:

A 1 = 1 0 0 0 µ1 cos(φ1) -µ1 sin(φ1) 0 µ1 sin(φ1) µ1 cos(φ1) , A 2 = µ2 cos(φ2) -µ2 sin(φ2) 0 µ2 sin(φ2) µ2 cos(φ2) 0 0 0 1 , (29)
where the numerical values of the parameters are µ 1 = 0.9, µ 2 = 0.2, φ 1 = π 6 , φ 2 = π 3 . Let us notice that ρ(A) = 1 with the Euclidean norm • being extremal for A.

We first compute a lower bound for the ρ-SJSR using Proposition 5. Let us remark that

N I = {A 1 A 2 , A 2 A 1 }.
We use the JSR toolbox [START_REF] Vankeerberghen | JSR: A toolbox to compute the joint spectral radius[END_REF] to compute tight bounds on ρ(N I ) and we obtain ρ(N I ) ∈ [ρ, ρ] where ρ = 0.7640321 and ρ = 0.7640322. It follows from Proposition 5, that

∀ρ > 1, λ(A, ρ) ≥ ρ ρ m .
We now apply Theorem 2 and Lemma 3 to compute an upper bound on the ρ-SJSR. Considering K = ρ, we search for an induced matrix norm

• * such that A 1 A 2 * ≤ ρ and A 2 A 1 * ≤ ρ.
Limiting the search to a quadratic norm of the form x * = x Qx, Q = Q ≥ 0 with the associated induced matrix norm M * = Q 1/2 M Q -1/2 , the conditions above are equivalent to the following LMIs: rem 3 translate to the following LMIs:

(A 1 A 2 ) QA 1 A 2 ≤ ρ 2 Q, (A 2 A 1 ) QA 2 A 1 ≤ ρ 2 Q.
I n ≤Q i , i ∈ I A i Q i A i ≤ρ 2 Q i , i, i ∈ I, i = i A i Q i+1 A i ≤ρ 2 Q i , i ∈ I \ {m} A m Q 1 A m ≤ρ 2 λ 2m Q m .
Similarly, searching for Lyapunov functions of the form V (J, x) = x Q J x, J I, the conditions of Theorem 4 translate to the following LMIs:

I n ≤Q J , J I A i Q J∪{i} A i ≤ρ 2 Q J , if J ∪ {i} = I A i Q 1 A i ≤ρ 2 λ 2 Q J , if J ∪ {i} = I.
Then, for each ρ ∈ Ω = {1, 1.05, • • • , 1.35}, we can find an upper bound λ(ρ) on λ(A, ρ) by searching (using a line search) for the smallest value of λ for which the LMIs above have a solution. Then, we can obtain an upper bound for the ρ-SJSR for all ρ ∈ (1, R) using ( 5):

∀ρ ∈ (1, R), λ(A, ρ) ≤ min ρ ∈Ω,ρ ≤ρ ρ ρ m λ(ρ ).
The resulting upper-bounds as well as the lower bound are shown in Figure 4. One can check that the upperbound computed using Theorem 4 is always tighter than that provided by Theorem 3. However, let us remark that the former approach involves solving a set of m(m + 1) LMIs while the latter requires solving a set of (2 m -1)(m+1) LMIs. For large values of m the latter approach is likely to be intractable.

Case study: synchronization over a network

We provide a short case-study regarding an application of our tools to the design of networked controllers for syn-chronization of distributed oscillators. Let us consider m + 1 identical components consisting of discrete-time oscillators whose dynamics is given by:

z i (t + 1) = Rz i (t) + u i (t), i = 1, . . . , m + 1 (30)
where

z i (t) ∈ R 2 , u i (t) ∈ R 2 and R = µ cos(φ) -µ sin(φ) µ sin(φ) µ cos(φ)
with µ = 1.02 and φ = π 6 . The input u i (t) is used for synchronization purpose and is based on the available information at time t. The information is exchanged over a single hop communication network, which consists of m communication channels between components i and i + 1, i = 1, . . . , m. We assume that the network capacity is limited, in such a way that only one channel can be activated at a given time instant. Formally, the active channel is given by a switching signal θ : N → I = {1, . . . , m}. Then, the input value implementing the synchronization protocol is given as follows:

u 1 (t) = k(z 2 (t) -z 1 (t)), if θ(t) = 1 0, otherwise u i (t) =    k(z i-1 (t) -z i (t)), if θ(t) = i -1 k(z i+1 (t) -z i (t)), if θ(t) = i 0 otherwise for i = 2, . . . , m, u m+1 (t) = k(z m (t) -z m+1 (t)), if θ(t) = m 0 otherwise (31)
where k ∈ (0, 1) is a control gain. Denoting the vector of synchronization errors as x(t) = (x 1 (t) , . . . , x m (t) ) with x i (t) = z i+1 (t) -z i (t), the error dynamics is described by a 2m-dimensional switched linear system of the form (1) with m modes. The expression of the matrices A i ∈ R 2m×2m , i = 1, . . . , m can be easily derived from ( 30)- [START_REF] Wang | Stability analysis of switched linear systems defined by regular languages[END_REF].

For m ≥ 2, it is clear that the synchronization cannot be achieved using arbitrary switching signals. Indeed, considering for instance the constant switching signal θ(t) = 1, for all t ∈ N, components 1 and 2 will synchronize but other agents that do not get the opportunity to exchange information cannot synchronize. Hence, we use shuffled switching signals. In our case study, we aim at providing answers to the following questions:

(1) Co-design a control gain k and a minimal shuffling rate γ * so that all oscillators synchronize; (2) Determine the maximal number of oscillators that can be synchronized given the proposed network architecture and capacity.

We compute ρ(A) using the JSR toolbox. We then use Theorem 4 to compute an upper bound λ(ρ) on the ρ-SJSR by solving the associated LMIs presented in the previous subsection. Note that

γ * := inf ρ>ρ(A) - ln(ρ) ln(λ(ρ)) ≥ inf ρ>ρ(A) - ln(ρ) ln(λ(A, ρ))
.

In this example, we observe numerically that the infimum is reached for ρ = ρ(A).

We computed γ * for several values of k ∈ (0, 1) and m + 1 ∈ {3, 4, 5} and we show the result in Figure 5.

In this figure, the dashed lines correspond to the maximal achievable shuffling rate 1 m . From Corollary 1, we know that the switched system can be stabilized if we use shuffled switching signals θ ∈ S s (I) whose shuffling rate γ θ ∈ (γ * , 1 m ]. Hence, in order to stabilize the system, we should carefully select the control gain k such that the corresponding value of γ * < 1 m . This answers the first point of our case study. For instance, we can see that for a system composed of m + 1 = 3 oscillators, for a control gain k = 0.6 we get γ * = 0.08 and we can stabilize the system by imposing a shuffling rate γ θ ∈ (0.08, 0.5]. If we choose k = 0.88, then γ * > 1 m and we are not able to stabilize the system by shuffling.

We also verified that for m + 1 = 6 oscillators, γ * > 1 m for all control gain k ∈ (0, 1). Therefore, for 6 oscillators, there is no suitable choice of control gain k such that the system can be stabilized by shuffling. This answers the second point of our case study.

We now proceed with some illustrative numerical simulations. Let us consider a system composed of m + 1 = 3 oscillators with control gain k = 0.4, then from Figure 5, the corresponding lower bound on the shuffling rate is γ * = 0.06. Let us consider the initial synchronization error x(0) = (-1.5 -0.5 2 -1) . We consider random switching signals generated by a discrete-time Markov chain with 2 states and the following transition matrix:

1-p p p 1-p where p ∈ (0, 1) is the probability to switch at a given time instant. It is easy to see that such switching signals are shuffled almost surely. Moreover, the larger p the higher the shuffling rate of the switching signals. We first consider p = 1/10, Figure 6 shows the evolution of κ θ (t) t , the switching signal θ(t) and the synchronization errors x 1 (t), x 2 (t). It is interesting to note that κ θ (t) t > 0.06 and that the system stabilizes as expected. Next, we consider p = 1/70, the corresponding simulation results are shown in Figure 7. We can check on the figure that in this case κ θ (t) t < 0.06 and that the system does not stabilize, the shuffling being too slow.

Conclusion

In this paper, we introduced the ρ-SJSR and the SJSR, a new notion of joint spectral radius for discrete-time switched linear systems driven by shuffled switching signals. We have established some of their properties and highlighted their relation to stability properties of switched systems. We presented a method based on the JSR of a finite set of matrices to compute asymptotically tight lower and upper bounds of the ρ-SJSR. We also developed two approaches based on Lyapunov functions and automata theoretic techniques to compute upper-bounds of the ρ-SJSR. While the first approach is computationally more tractable, the second approach was shown to be tight. Numerical examples show the effectiveness of the proposed numerical approximation methods. We also showed an application of the SJSR to the synchronization of unstable oscillators.

The current work opens several research directions for the future. Firstly, even though the proposed Lyapunov conditions make it possible to compute tight upper bounds of the ρ-SJSR, their translation to linear matrix inequalities may introduce some conservatism. It is then interesting to investigate how techniques similar to the path-complete graph Lyapunov functions [START_REF] Ahmadi | Joint spectral radius and path-complete graph lyapunov functions[END_REF] can be used in order to derive linear matrix inequalities whose solution provide a tight upper bound of the ρ-SJSR. Secondly, it has been shown in [START_REF] Wang | Stability analysis of switched linear systems defined by regular languages[END_REF] that stability of a switched system with switching signals belonging to an ω-regular language is equivalent to the shuffled stability of a lifted switched system, so our current approach can readily be used to analyze stability properties of such systems. However, since our approach is based on automata theoretic techniques, it is natural to think that one can derive stability conditions by working directly on the Büchi, Rabin or Muller automaton specifying the ω-regular language.

Fig. 2 .

 2 Fig. 2. Automaton describing the dynamics of η in (24) for m = 3. State labels correspond to the value of η, transition labels correspond to the value of θ.

Fig. 3 .

 3 Fig. 3. Automaton describing the dynamics of η in (28) for m = 3. State labels correspond to the value of η, transition labels correspond to the value of θ.

.Fig. 4 .

 4 Fig. 4. Lower and upper bounds on the ρ-SJSR for (29).

Fig. 5 .

 5 Fig. 5. Minimal shuffling rate γ * as a function of the control gain k and of the number of oscillators m + 1.

  Fig. 6. Time evolution of the synchronization error x(t)the switching signal θ(t) for a shuffling rate higher than γ * = 0.06. Fig.7. Time evolution of the synchronization error x(t) and the switching signal θ(t) for a shuffling rate lower than γ

	(t)/t	0.1 0.2
				0
					0	100	200	300	400	500
				2
	(t)	1.5
				1
					0	100	200	300	400	500
	x 12 (t) 11 (t), x	-2 0 2	0	100	200	300	400	500 x 11 (t) x 12 (t)
	x 22 (t) 21 (t), x	-2 0 2	0	100	200	300	400	500 x 21 (t) x 22 (t)
							t
				0.1
	(t)/t	0.05
					0
					0	100	200	300	400	500
					2
	(t)	1.5
					1
					0	100	200	300	400	500
	x 12 (t) 11 (t), x	-400 -200 0 200 400	x 11 (t) x 12 (t)
					0	100	200	300	400	500
	x 22 (t) 21 (t), x	-200 0 200	x 21 (t) x 22 (t)
					0	100	200	300	400	500
							t

* = 0.06.

The Matlab scripts of the two numerical examples are available at the following repository: https://github.com/ georgesaazan/Shuffled-systems

Georges Aazan. This project has received funding from the H2020-EU.1.1. research and innovation programme(s) -ERC-2016-COG under grant agreement No 725144. This work was supported in part by the Agence Nationale de la Recherche (ANR) under Grant HANDY ANR-18-CE40-0010.


Guillaume Aoust
email: guillaume.aoust@onera.fr

Raphaël Levy

Béatrice Bourgeteau

Olivier Le Traon

Acoustic damping on flexural mechanical resonators

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Flexural mechanical resonators are most of the time operated in vacuum in order to minimize any coupling with the surroundings that could affect their excellent resonance characteristics. Some applications however require a complete immersion within a fluid medium. For example, they are used for highly sensitive chemical sensing platform [START_REF] Li | Nanoelectromechanical resonator arrays for ultrafast, gas-phase chromatographic chemical analysis[END_REF], but also for in-situ force measurements such as the photo acoustic detection [START_REF] Kosterev | Quartz-enhanced photoacoustic spectroscopy[END_REF], the resonant optothermoacoustic detection [START_REF] Kosterev | Resonant optothermoacoustic detection: technique for measuring weak optical absorption by gases and micro-objects[END_REF] or precise temperature measurements [START_REF] Hopcroft | Using the temperature dependence of resonator quality factor as a thermometer[END_REF].

In those cases, the surrounding fluid strongly modifies the characteristics of the resonance and can dramatically reduce the sensors performances compared with vacuum operation. It is now established that the damping occurring within the fluid can be attributed to several mechanisms, among which the most relevant ones are viscous damping and acoustic radiation damping. These two effects, originating both from a fluid-structure interaction, rely on two distinct physical properties of the fluid that are viscosity and compressibility, and their contributions behave differently with respect to the dimensions of the resonator.

Acoustic radiation damping is not frequently considered in high quality flexural resonators. One of the main reasons is that a huge majority of them are designed for applications in the inertial or time-frequency domains, i.e. under vacuum operation where no fluid-structure occurs. For other applications requiring an immersion within a fluid, viscous damping has been shown to be negligible for macroscopic resonators because of their low surface over volume ratios. Their acoustic properties have been investigated to obtain reliable expressions for the total quality factor [START_REF] Blake | On the damping of transverse motion of free-free beams in dense stagnant fluids[END_REF][START_REF] Jeyapalan | Radiation efficiencies of beams in flexural vibration[END_REF][START_REF] Johnston | Acoustic and internal damping in uniform beams[END_REF]. For smaller resonators, the number of geometries reported up to now remains quite limited; and acoustic damping is most of the time hardly noticeable if not negligible. Nevertheless, acoustic damping is still possible with specific geometries and can become comparable with viscous damping; therefore a more detailed study is of significant importance for any optimization purposes.

Finite element simulation softwares can quantitatively and accurately predict the mechanical behavior of any type of resonator. However, we will focus on analytical models in order to keep a good physical insight. In the following, we will restrain to the case of resonators made of beams with a circular or rectangular cross section, since it almost covers all kind of existing structures.

A first complete analytical model including both viscosity and compressibility has already been developed [START_REF] Van Eysden | Small amplitude oscillations of a flexible thin blade in a viscous fluid: exact analytical solution[END_REF]. The result is exact but only deals with infinitely thin cantilevers and fails to provide either a direct physical insight of the compressibility effect or a simplified formula for the quality factor. Moreover, the approach does not include a possible static wall nearby, which is necessary to deal with viscous damping for resonators made of several beams such as tuning forks [START_REF] Aoust | Viscous damping on flexural mechanical resonators[END_REF].

Other existing analytical treatments currently available in the literature generally assume that the surrounding fluid acts as an uncompressible medium [START_REF] Aoust | Viscous damping on flexural mechanical resonators[END_REF][START_REF] Cox | Characteristics of laterally vibrating resonant microcantilevers in viscous liquid media[END_REF]. The damping is then reduced to the fluid viscous effects, for which analytical results give quite accu- rate formulas for the quality factor. However, such an assumption is not valid in some cases. It has been evidenced in the study of tuning forks in gas [START_REF] Christen | Air and gas damping of quartz tuning forks[END_REF], and more recently in cryogenic fluids [START_REF] Tuoriniemi | Mode analysis for a quartz tuning fork coupled to acoustic resonances of fluid in a cylindrical cavity[END_REF][START_REF] Bradley | Crossover from hydrodynamic to acoustic drag on quartz tuning forks in normal and superfluid He[END_REF]]. In the latter example, it has been shown that commercial quartz tuning forks immersed within liquid or gaseous helium at very low temperature can provide various parameters such as viscosity, pressure or temperature [START_REF] Blaauwgeers | Quartz tuning fork: thermometer, pressure-and viscometer for helium liquids[END_REF]. Under such conditions, acoustic damping is not a negligible and has to be taken into account.

From a theoretical point of view, both viscous and acoustic damping effects are accounted for in the Navier-Stokes equations. However, no general analytical solution has yet appeared if both sources of damping are considered simultaneously. Indeed, the coupled equation system governing the behavior of the fluid on the first hand and the motion of the resonator on the other hand cannot be reduced into two independent equations.

A simplified approach consists in solving the two problems in a sequential way, which is an approximate way to force the independence. The motion of the resonator is solved first assuming that the surrounding fluid is uncompressible, leading to an analytical expression of the damping only due to viscous effects. Then, the motion of the resonator is used as a boundary initial condition for the purely acoustic problem neglecting viscosity.

In this paper, we will generalize this sequential approach to any flexural mechanical resonator made of beams with arbitrary rectangular or circular cross section, as well as investigate the validity of our results with numerical simulations and experiments. We also provide expressions for the quality factors, as well as simplified expressions to capture more easily the physics behind fluid damping on resonators.

Theoretical model

The single beam model

We first consider a single beam oscillating within a homogeneous fluid unbounded in space. Notations used in the following are presented in Fig. 1. The motion of the resonator subject to viscous damping only has already been solved using the Euler-Bernoulli bending theory, so we will use directly the results obtained in this previous work [START_REF] Aoust | Viscous damping on flexural mechanical resonators[END_REF]. Let us remind the main assumptions concerning the beam, which generally hold true for any flexural resonator:

• The beam's cross section is rectangular and uniform over its entire length L. • The length of the beam L greatly exceeds its transverse dimensions e and l. • We consider only flexural modes of vibration in the y-direction, whose amplitudes are supposed far smaller than dimensions e and l. Plane sections hence remain plane and normal to the axis of the beam. • The material composing the beam is assumed isotropic and homogeneous. Its young modulus E and density b are constants.

Considering the surrounding fluid, we suppose that:

• The fluid is incompressible with a homogeneous density f and dynamic viscosity . • The fluid is in the continuum regime [START_REF] Li | Review of viscous damping in micro-machined structures[END_REF] and unbounded in space.

The latter regime is defined by its Knudsen number K n smaller than 0, 01, ensuring that the Navier-Stokes equation are valid.

Under the previous assumptions, the position of the oscillating beam Ŵ (x, ) in the y direction can be written on its n th eigen angular frequency ω n as Eq. [START_REF] Li | Nanoelectromechanical resonator arrays for ultrafast, gas-phase chromatographic chemical analysis[END_REF].

Ŵ = Ŵ (ω) n (x) = L 0 dx F drive (x, ω) n (x) 12jL 4 Q v (ωn) Ee 3 l˛4 n L 0 2 n (x) dx n (x) (1)
where ˛n and n are respectively the characteristic coefficient and the normalized deformation of the beam when oscillating on its n th mode of vibration. Their values can be found in any textbook about dynamic beam flexure, and two useful cases are recalled in Table 1.

The viscous quality factor Q v has been derived in the previous study [START_REF] Aoust | Viscous damping on flexural mechanical resonators[END_REF]. Function F drive is the excitation force acting on the beam in the y direction other than the viscous forces.

In this second step, we suppose that the motion of the beam expressed in Eq. ( 1) is given. The strategy to derive the acoustic

˛n = n + 3 2 (n > 2)
quality factor Q a is to solve first the wave equation governing the pressure field all around the resonator so as to obtain the acoustic radiated power P.

Introducing the wavenumber k = ω/c, the pressure field p around the resonator in a 3D space is the solution of problem (P)

(P) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∇ 2 p + k 2 p = 0 → ∇ • → n p = f ω 2 Ŵ → y • → n lim r→+∞ r ∂ ∂r -jk p = 0 (2)
The first Eq. in (2) is the classic Helmholtz equation where any viscous damping is neglected, whereas the second equation is the boundary conditions at the fluid-resonator interface of normal vector → n . The third equation is known as the Sommerfeld radiation condition and set the boundary condition far from the resonator. The boundary conditions on the resonator will impact the acoustic quality factor depending on the number of beams composing the resonator. For a single beam oscillating in an unbounded fluid medium, two different kinds of analytical treatments have been proposed in the past.

The acoustic equivalent approach

The simplest approach establishes an analogy between the oscillating beam and a set of acoustic point sources, as depicted in Fig. 1 (b). Since any beam's cross section communicates a velocity to the fluid in the direction of its vibration, it therefore constitutes an acoustic source. It is known that the pressure field due to a point source located at r = r 0 is solution to the following acoustic problem [START_REF] Morse | Methods of Theoretical Physics[END_REF]:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∇ 2 p + k 2 p = j f ωq (r 0 ) ␦ (r -r 0 ) q (r 0 ) = → U • → n dS lim r→+∞ r ∂ ∂r -jk p = 0 (3)
In the previous equation, q is the volume flow rate, defined within the analogy by the integral on the surface close to r 0 of the beam's normal velocity

→ U • → n .
The total contribution of the different beam's sections to the acoustic radiated power is however not straightforward. Two ways of including this dependency can be used.

2.1.1.1. The point source approach. The first way is to consider two point sources in a 3D space, and we will refer to this technique as the point source approach. The system is equivalent to a single section of the beam carrying the integrated acoustic strength over the whole beam. The volume flow rate q is set as the total volume flow rate integrated along the beam, and will take at least roughly into account the possible interferences in the pressure field due to the length of the beam. To do so, we will add a dephasing depending on the polar angle ϕ. An example of this method can be found in the work of Jünger, who however restricted his approach to the case ke << 1 [START_REF] Junger | Sound radiation by resonances of free-free beams[END_REF]. We hence choose the volume flow rate displayed in Eq. ( 4) to make the analogy.

q = jωl Ŵ (ω) L 0 n( x) e -jkxcosϕ dx (4)
We use the spherical coordinates r, ϕ, pictured in Fig. 1 to solve the problem. The 3D Green's function associated with the wave equation on pressure being G (r, r 0 , ω) = exp (jk|r -r 0 |) / (4|r -r 0 |), we obtain the solution for the pressure by direct convolution:

p = -jω f q 4 e jkR+ R + - e jkR- R - (5)
The numbers

R ± = r 2 + e/2 2
± re cos measure the distances between the current point of observation and the two point sources. In order to separate the coordinates, we use the deported potential identities valid if r > e/2 [START_REF] Morse | Methods of Theoretical Physics[END_REF]:

e jkR± R ± = jk +∞ m=0 (2m + 1) P m cos j m ke 2 h m (kr) (6)
In the previous expression, P m are the Legendre polynomials, j m and h m are respectively the spherical Bessel and Hankel functions.

Gathering the two source terms within the same infinite sum, we notice that the difference P m cos -P m -cos appears and cancels out for any even integer m because of the parity properties of the Legendre polynomials. Eq. ( 5) can hence be rewritten as Eq. ( 7).

p = -jω f q 4 2jk +∞ m=0 (4m + 3) P 2m+1 (cos )j 2m+1 ke 2 h 2m+1 (kr) (7)
The power radiated by the resonator far from the resonator is given by Eq. [START_REF] Van Eysden | Small amplitude oscillations of a flexible thin blade in a viscous fluid: exact analytical solution[END_REF].

P (ω) = =0 2 ϕ=0 |p R, ϕ, | 2 2 f c r 2 sin ddϕ (8)
Using the ortho-normality property of the Legendre Polynomials 0 P i cos P j cos sin d = 2ı ij / (2i + 1) (where ı ij represents the Kronecker symbol), and taking the asymptotic limit of the spherical Hankel functions far from the resonator, we hence obtain from Eqs. ( 7) and ( 8) the expression of the radiated acoustic power.

P (ω) = ω 4 l 2 f 4 2 c Ŵ (ω) 2 2 ϕ=0 L 0 n( x) e -jks cos ϕ dx 2 dϕ +∝ m=0 (4m + 3) j 2m+1 ke 2 2 (9)
If E (ω) denotes the energy stored in the resonator, then the quality factor can be defined as:

Q a = ωE (ω) P (ω) (10)
If we suppose that the energy stored in the beam is only kinetic energy, we have according to the Euler-Bernoulli formalism

E (ω) = 1 2 Ŵ (ω) 2 b elω 2 L 0 2 n (x) dx (11)
Combining Eqs. ( 9), ( 10) and [START_REF] Christen | Air and gas damping of quartz tuning forks[END_REF], we obtain the acoustic quality factor of the beam:

Q a = 2 2 b f e kl 0 L 2 n (x) dx ϕ=0 2 0 L n (x) e -jkx cos dx 2 d +∞ m=0 (4m + 3) j 2m+1 ke 2 2 (12)
2.1.1.2. The infinite line source approach. The second approach has been proposed by Schmoranzer [START_REF] Schmoranzer | Acoustic emission by quartz tuning forks and other oscillating structures in cryogenic He fluids[END_REF] in the case of a clamped free tuning fork. His approach can be extended to the single beam with any boundary conditions as well. We will refer to this method as the infinite line source approach. His result concerning tuning forks will be reported in the next section.

The idea is to introduce from the start the interferences caused by the different sections of the beam in a more rigorous manner. In order to do so, the beam is supposed infinite in the x direction, and the acoustic equivalent problem is now represented by two infinite line sources of opposite strengths. The resulting translation invariance simplifies the problem (P) to a 2D problem, and the volume flow rate is now given per unit length in Eq. [START_REF] Bradley | Crossover from hydrodynamic to acoustic drag on quartz tuning forks in normal and superfluid He[END_REF].

q = jωl Ŵ (ω) (13)
Using the cylindrical coordinates r, , x (which is equivalent to the previous spherical coordinates r, 0, in the y O z plane, taking x as the Cartesian coordinate along the x axis) the Green function is now in 2D and reads as G (r, r 0 , ω) = jH 0 (k|r -r 0 |) /4, with H 0 the Hankel function of the first kind. All of the previous steps are then followed and adapted to this new situation, therefore we will just give here the final expression. The new pressure field is

p = ω f q 4 [H 0 (kR+) -H 0 (kR-)] (14)
The new deported potential identities valid if r > e/2 [START_REF] Morse | Methods of Theoretical Physics[END_REF] are now:

H 0 (kR) = m=+∞ m=-∞ e jm(-0) H m (kr) J m (kr 0 ) ( 15
)
where H m and J m are now respectively the Hankel and Bessel functions of the first kind, and r 0 , 0 designates the coordinates of the center of the line source. The total pressure hence reads:

p = ω f q 2 m=+∞ m=-∞ e j(2m+1) H 2m+1 (kr) J 2m+1 ke 2 ( 16
)
The approximation of this method comes during the integration of the power radiated. To obtain an analytical expression, any power radiated in the half spaces above and below the beam's end is neglected. The integration can then only occur far from the resonator where it is possible to use the asymptotic expression for the Hankel function. Also, integration is needed along the length of the beam to account for the finite dimensions of the resonator. We recover roughly the dependency of the volume flow rate using q = jωl Ŵ (ω) n (x), which is a new approximation because we had already assumed invariance along the x axis.

We obtain the radiated power:

P (ω) = ω 3 f Ŵ (ω) 2 l 2 L 0 2 n (x) dx m=+∞ m=0 J 2m+1 ke 2 2 ( 17
)
We deduce immediately the corresponding quality factor.

Q a = 1 2 b f e l 1 m=+∞ m=0 J 2m+1 ke 2 2 (18)
2.1.1.3. Conclusion. None of the two previous models is completely rigorous concerning the role of the beam's length on the quality factor. Both try to capture its influence on the radiated power, and a comparison is needed to evaluate which one shows the more accurate results. We expect however that the line source method will be better if the beam's length is large compared to the acoustic wavelength (L >> ) since it matches well the initial approximation of invariance along the length. Similarly, we expect the point source method to be better if the beam's length is short compared to the wavelength (L << ), since it reduces the impact of the dephasing along the beam that was added in a non-rigorous manner.

The beam cross section approximation

The second kind of formalism makes an assumption on the shape of the beam's section in order to obtain an analytical solution of the original vibroacoustic problem. This solution seems more rigorous since no analogy is used, and has been proposed by modeling the beam's section as a circular section beam [START_REF] Johnston | Acoustic and internal damping in uniform beams[END_REF] or an ellipsoidal section beam [START_REF] Johnston | Acoustic and internal damping in uniform beams[END_REF][START_REF] Morse | Methods of Theoretical Physics[END_REF][START_REF] Blake | The radiation from free-free beams in air and in water[END_REF]. However, additional assumptions are required.

To model arbitrary rectangular cross sections, the ellipsoidal approach seems closer than the cylinder in terms of shape. In the elliptic cross section modeling, the use of elliptic coordinates allows to derive a solution written in terms of Mathieu functions. We will refer to this method as the elliptic approach. The formalism however requires the beam's length L to be much greater than the acoustic wavelength = 2c/, because the assumption of an infinitely long beam is needed. If those models are used on the lowest flexural mode number, discrepancies are expected since the edge effects become important.

Since the total radiated power has been derived several times in the literature, we will give directly the expression given by Blake [START_REF] Blake | The radiation from free-free beams in air and in water[END_REF] assuming ke << 1:

P (ω) = f c 512 (kl) 4 ω 2 ϕ=0 sin 3 ϕ L 0 n( x) e -jkx cos ϕ dx 2 dϕ Ŵ (ω) 2 (19)
Using the energy stored in the resonator E (ω), we directly obtain the quality factor:

Q a = 256 b f 1 (kl) 3 × e 0 L 2 n (x) dx ϕ=0 sin 3 ϕ| 0 L n (x) e -jkx cos ϕ dx| 2 dϕ ( 20
)
The result is similar to the expression given by Lochon et al. [START_REF] Lochon | A microcantilever chemical sensors optimization by taking into account losses[END_REF] divided by a 4 2 factor they omitted.

Total quality factor

We write the total quality factor of the beam vibrating in the compressible medium as the sum of its different contributing effects we just obtained sequentially, according to the usual additive rule valid for high quality factor resonators:

1 Q = 1 Q v + 1 Q vac + 1 Q a (21)
• Q vac is the quality factor of the beam in the absence of the surrounding fluid. This value will depend on the anchor and internal losses, and are typically much lower than those due to the interaction with the fluid. • Q v is the quality factor of the beam due to the viscosity of the surrounding fluid only, considering that the vacuum quality factor is infinite and that the fluid is incompressible. We obtained its expression in [START_REF] Aoust | Viscous damping on flexural mechanical resonators[END_REF]. • Q a is the quality factor associated with the acoustic losses, assuming that the vacuum quality factor is infinite and that the fluid is inviscid. We obtained its expression in the previous sections.

All the presented formulas of acoustic losses will be used and compared in the next section of this work.

The tuning fork case

We consider now a tuning fork immerged in a compressible unbounded fluid. We will apply the same strategy to obtain the expression of the total quality factor.

It is known that the viscous damping in an incompressible medium is modified for the tuning fork by adding a new contribution Q tf , accounting for [START_REF] Aoust | Viscous damping on flexural mechanical resonators[END_REF]:

• The squeeze film damping appearing because of the proximity of the prongs if the tuning fork is vibrating in an in-plane mode of vibration. • The Couette flow damping appearing because of the proximity of the prongs if the tuning fork is vibrating in an out-of-plane mode of vibration.

In any case, the total viscous damping quality factor can be written as the sum of the tuning fork specific losses (squeeze film damping and Couette flow damping) and of the lateral and front viscous losses as expressed in [START_REF] Aoust | Viscous damping on flexural mechanical resonators[END_REF], and reads

1 Q v = 1 Q tf + 1 Q lat + 1 Q front ( 22
)
Results obtained in the single beam case cannot be adapted to the tuning fork with acoustic damping. However, the same method can be used to derive the new results.

The point source approach and the infinite line approach for inplane tuning fork vibration have been first obtained respectively by Sillitto [START_REF] Sillitto | Angular distribution of the acoustic radiation from a tuning fork[END_REF] and Schmoranzer et al. [START_REF] Schmoranzer | Acoustic emission by quartz tuning forks and other oscillating structures in cryogenic He fluids[END_REF]. We will not deal with the analogy of the infinite cylinders [START_REF] Clubb | Quartz tuning fork viscometers for helium liquids[END_REF] since this approach is very similar to that of Schmoranzer.

In-plane motion: Sillitto extended model

We will use the same notations as in the single beam case, except that the origin is relocated between the prongs as showed in Fig. 2. Since the tuning fork is composed of two beams, the acoustic analogy becomes a set of four acoustic point sources. The model has been first proposed by Sillitto only for an in-plane mode of vibration, and we will slightly extend it in an attempt to overcome the L << assumption he uses by adding a rough interference contribution.

For an in-plane mode of vibration, the pressure is now given by

p = -ω 2 f qe -jωt 2c +∞ m=0 (4m + 1) P 2m (cos ) j 2m k g 2 + e -j 2m kg 2 h 2m (kr) (23)
Using Eq. ( 8), the total radiated power far from the resonator is:

P (ω) = f ω 2 2c |q| 2 +∞ m=0 (4m + 1) j 2m k g 2 + e -j 2m kg 2
Table 2 Simplified expressions for the acoustic quality factor obtained using different approaches. The validity conditions of the formulas are detailed in the column "Hypothesis".

Approach Hypothesis Expression of Qa

Single beam

Point source

ke 2 ≪ 1 24 2 b f 1 k 3 el 0 L 2 n (x)dx ϕ=0 2 | 0 L n (x)e -jkx cos ϕ dx| 2 dϕ ke 2 ≪ 1 L ≪ 12 b f 1 k 3 el 0 L 2 n (x)dx 0 L n (x)dx 2 Infinite line ke 2 ≪ 1 8 b f 1 k 2 el Elliptic ke 2 ≪ 1 256 b f 1 (kl) 3 e 0 L 2 n (x)dx ϕ=0 sin 3 ϕ| 0 L n (x)e -jkx cos ϕ dx| 2 dϕ ke 2 ≪ 1 L ≪ 192 b f 1 (kl) 3 e 0 L 2 n (x)dx 0 L n (x)dx 2
Tuning fork

Inplane

Point source

k e + g 2 ≪ 1 80 2 b f 1 k 5 el[e + g] 2 0 L 2 n (x) dx ϕ=0 2 | 0 L n (x) e -jkx cos ϕ dx| 2 dϕ k e + g 2 ≪ 1 L ≪ 40 b f 1 k 5 el(g+e) 2 0 L 2 n (x)dx 0 L n (x)dx 2 Infinite line k e + g 2 ≪ 1 64 3 b f 1 k 4 el(e+g) 2 Out-of-plane Infinite line k 2 (l + g) 2 + e 2 ≪ 1 64 b f 1 k 4 el(l+g) 2
We choose to use the same expression as we did in the single beam case, defined in Eq. ( 4).

By noticing that the energy stored in the tuning fork is now twice that contained in one of its prongs whose expression has already been given in Eq. ( 11), we obtain the following quality factor: We retrieve the case developed by Sillitto by assuming that the total length is small compared to the acoustic wavelength: the dephasing is then negligible and the factor B ≡ q used by Sillitto becomes B = jl Ŵ (ω) L e , with L e = L 0 n (x) dx the effective length of the prongs.

Qa = 4 2 b f e kl 0 L 2 n (x) dx

In-plane motion: Schmoranzer model

A different approach to overcome the L << has been recently proposed by Schmoranzer. We already exposed his technique in the single beam case, we will hence here only give the main expressions. The power radiated and integrated on a cylinder far from the sources, neglecting the half spaces above and below the resonator, is:

P (ω) = f ω 3 Ŵ (ω) 2 l 2 L 0 2 n (x) dx 2 m=+∞ m=-∞ J 2m k e + g 2 -J 2m kg 2 2 (26)
And the corresponding quality factor using the Kronecker notation ı 0,m is:

Q a = 2 b f e l 1 m=+∞ m=0 2 1+ı 0,m J 2m k e + g 2 -J 2m kg 2 2 (27)

Out-of-Plane motion: Schmoranzer extension

Finally and for theoretical purposes only, we provide an extension of the model proposed by Schmoranzer for the tuning fork on an out-of-plane mode of vibration. We draw the reader's attention to the fact that the dimensions of the fork e and l are switched in this case, in order to keep the formulas for viscous damping unchanged (see Fig. 2(b). The power radiated and integrated on a cylinder far from the sources, neglecting the half spaces above and below the resonator, reads:

P (ω) = 4 3 f Ŵ (ω) 2 l 2 L 0 2 n (x) dx m=+∞ m=1 sin(2m 0 )J 2m (kr 0 ) 2 (28)
With 0 = tan -1 e/ (g + l) and r 0 = (l + g) 2 + e 2 /2. The corresponding quality factor is:

Q a = 1 4 b f e l 1 m=+∞ m=1 sin 2m 0 J 2m (kr 0 ) 2 (29)

Total quality factor

The same expression of the total quality factor we obtained for the single beam stands for the tuning fork:

1 Q = 1 Q v + 1 Q vac + 1 Q a ( 30
)
Tuning forks typically present excellent decoupling from the support, minimizing the support losses and ensuring a high vacuum quality factor, often greater than Q vac ∼10 5 .

Summary of the theoretical part

We also provide simplified expressions for the previous acoustic quality factors, valid when the acoustic wavelength is large compared to the width of the resonators, and obtained using Taylor expansion over the infinite summations (Table 2).

We will now investigate the validity of the approach combining both viscous and acoustic damping obtained independently, with experiments where neither the viscous nor the acoustic dampings are negligible. We will show that finite element simulation softwares can give a precise value for the radiation losses and in very good agreement with the experiments. We will then compare the results obtained with our various acoustic analytical models.

Experiments

All the following results were obtained considering the fluid medium as air at room temperature. In order to investigate their pressure dependency, the resonators are inserted inside a sealed vacuum bell. A pressure sensor TPR280 has been used to monitor the pressure variations. The vacuum bell is pumped through a sluice gate, and any pressure in the range1Pa -10 5 Pa can be achieved. The vacuum quality factor Q vac used in the modeling is considered to be the measured quality factor of the resonator when the pressure is lower than0.1Pa. The quality factors at various pressures are extracted from the measurement of the 3 dB band pass of the resonators piezoelectric response. We recall that the models developed in the previous sections assume that the surrounding fluid is characterized by a Knudsen number lower than 0.01. We have reported in the following figures the successive Knudsen regimes, indicating the pressure ranges where the models are supposed to be valid.

Numerical simulations have been performed with a finite element model developed using OOFELIE::Multiphysics © , where non-viscous elements have been chosen to simulate an inviscid and hence purely vibro-acoustic problem [START_REF] Paquay | Développement d ' une méthodologie de simulation numérique pour les problèmes vibro-acoustiques couplés intérieurs/extérieurs de grande taille[END_REF][START_REF]OOFELIE::Multiphysics©, User Manual[END_REF]. The solver type is harmonic and we used a perfectly matched layer (PML) as a boundary condition around the resonator to simulate acoustic waves propagating toward infinity. The numeric acoustic quality factors are hence defined by the ratio Q a = / , with ω the 3 dB band pass.

Single beam

We use as a first experiment a clamped-clamped quartz beam to validate the combined approach. The beam's features are l = 0.6mm, e = 0.2mm, L = 4mm and is operated on its first flexural mode of vibration at a frequency close to f = 69kHz. We obtained the results presented in Fig. 3.

It has been shown on a wide range of resonators that viscous damping can be predicted within a ±15% precision interval [START_REF] Aoust | Viscous damping on flexural mechanical resonators[END_REF]. Considering viscous damping only (orange dashed curve), the experiments displayed in Fig. 3(a) show a prediction error of approximately90%, and constitutes a strong indication that an additional damping is present. We used both viscous and acoustic analytical models developed in the previous section to predict the final quality factor of the resonator. We see that the combination of an analytical prediction of the viscous damping and a numerical prediction of the acoustic damping using the finite elements method is in very good agreement with the experiment. The predicted result discrepancies do not exceed6%, which is already below the prediction accuracy of the viscous damping analytical model. Considering this geometry of clamped-clamped beams, it seems that finite elements numerical simulations predict the acoustic damping quite accurately, and can therefore be used as a reference for accuracy comparisons of analytical acoustic damping models.

The fully analytical combination (the green plain line) has been obtained using the elliptic approach as the analytical model for acoustics. It leads to a total quality factor prediction within 35% predictive error. This value is impacted by the accuracy of the viscous damping model, the accuracy of the acoustic damping model as well as the solving technique implying to compute the two dampings mechanisms sequentially.

In Fig. 3(b), we see that all analytical predictions of the acoustic damping are under-evaluating the actual value, which should be somewhere close to the results obtained from the finite element simulations (black squares). However, each model displays the correct dependency with respect to the pressure, evolving with the density of the fluid as 1/ f . The two models involving an acoustic equivalence are largely under predicting the correct value by a factor of 10 for the infinite line source approach and 20 for the punctual source approach. On the contrary, the elliptic approach gives almost satisfying predictions with a factor 3.

Full numerical simulations taking into account both viscous and acoustic dampings would probably help to clarify the accuracy of the combination method. In any case, it seems that analytical formulas are not very accurate to predict the quality factor of single beams, at least on their lowest order of vibration as it is the case in our experiment.

Tuning forks

In-house tuning fork

We fabricated a custom tuning fork, which seems also subject to acoustic damping. We carefully designed the anchors, ensuring reduced anchor losses and providing an excellent vacuum quality factor of Q vac = 10 6 . The tuning fork is clamped-free and vibrates at approximately f = 42kHz on its first in-plane mode of vibration, and the tuning fork thickness l = 1mm has been chosen relatively large compared to the acoustic wavelength = 8mm in order to favor acoustic damping. The total length L = 6.8mm is of the order of the wavelength, and we have a prong width of e = 4mm.The results are displayed in Fig. 4.

Using the same argument as in the single beam case, the experimental data in Fig. 4(a) shows dramatic damping that cannot be explained by viscous effects only. Under the hypothesis that acoustic damping is responsible, we obtain very accurate prediction results for the method using finite elements approach. In this case, the analytical models also predict particularly well the damping. In Fig. 4(b), we see that the results obtained in Section 2.2.1 for the point source approach and in 2.2.2 for infinite line source approach are extremely close.

This experiment, which is performed in the case L∼, shows that the two analytical models are accurate and give almost the same results as the numerical simulations.

The commercial QTF discussion

It is also possible to apply the previous formulas to different tuning forks. A large majority of them used in the literature are commercial quartz tuning forks which are generally assumed to be only subject to viscous damping. Their typical dimensions are l = 0.34mm, e = 0.6mm, L = 3.75mm and g = 0.31mm, and they are most of the time operated on their first flexural mode of vibration at a frequency of approximately f = 32kHz.

The finite element software computes that, at atmospheric pressure in air, such a tuning fork displays an acoustic quality factor of Q a = 3.7 10 6 . This is much greater than the typical viscous quality factors 10 4 , which is why it is generally and legitimately assumed that the acoustic damping is negligible.

For those tuning forks, we have ke << 1 and it is possible to use the simplified expressions for the acoustic quality factors when it is operated on their first in-plane mode of vibration. The ratio between the point and the infinite line approach is:

Q point Q line ∼ 10 kL (31)
This simplified ratio shows that no unique analytical model has appeared to predict accurately the quality factor for any kind of shape, and the two methods scale differently with the fre- quency and the prongs length. The numerical applications give Q a = 2.3 10 6 and Q a = 1.1 10 7 respectively for the infinite line source approach and the point source approach for the standard tuning fork. From our experiment and given the assumptions made to derive those models, we expected the point source approach to be more accurate when the prong length is much shorter than the wavelength L << ; and the infinite line approach will be better otherwise. In this case, ∼10.6mm > L but the infinite line approach seems to produce a better prediction. When the acoustic wavelength and the prongs length are not too far apart, no conclusions can be drawn about picking one model over the other.

Conclusion

We showed that a sequential treatment of losses can accurately predict the quality factor of any flexural resonator in a fluid medium. The use of point sources as an acoustic equivalent for the tuning forks produces a very good value for the acoustic damping when the length of the prongs is roughly equal to the sound acoustic wavelength it generates. When the acoustic wavelength is greater than the prongs length, one should prefer the infinite line source approach. Concerning the single beam acoustic losses predictions, we showed that acoustic equivalent models can dramatically fail, whereas the elliptic approach gives a reasonable estimate. In any case and if the studied geometry does not fit the analytical model we described, it is still preferable to use finite element simulations software using non viscous acoustic medium to obtain a precise answer. Those results can be employed for the design of an optimized geometry, or example to obtain the best possibly achievable quality factor in fluid.

Fig. 1 .

 1 Fig. 1. Schematic of the acoustic equivalence for the single beam resonator. (a) The geometrical cross section of the oscillating beam of width e and thickness l. The amplitude of the displacement along the y axis is denoted Ŵ (x) . (b) The corresponding point sources modeling of the beam, with S the strength of the source.

Fig. 2 .

 2 Fig. 2. Schematic of a tuning fork and its acoustic point sources equivalence. (a) The tuning fork operated on an in-plane mode of vibration, and (b) the point source equivalent acoustic model. (c) The tuning fork operated on an out-of-plane mode of vibration and (d) the point source equivalent acoustic model.

Fig. 3 .

 3 Fig. 3. (a) Experimental investigation of the total quality factor of a single beam resonator subject to both viscous and acoustic damping. The analytical acoustic model used is the elliptic approach. (b) Comparisons between the different analytical acoustic damping approaches and finite element simulations. The analytical values of the viscous damping are also provided to compare the relative amplitudes.

Fig. 4 .

 4 Fig. 4. (a) Experimental investigation of the total quality factor of our homemade tuning fork subject to both viscous and acoustic damping. The analytical acoustic model used is the point source approach. (b) Comparisons between the different analytical acoustic damping approaches and finite element simulations. The analytical values of the viscous damping are also provided to compare the relative amplitudes.

Table 1

 1 Characteristic coefficient ˛n and deformation n of a beam vibrating on its n th flexural mode of vibration.

	Boundary conditions	Approximated Values	Deformation n
	Clamped-free	˛0 = 1.875		sin ˛n x L -sinh ˛n x L + cos ˛n +cosh ˛n sin ˛n -sinh ˛n	cos ˛n x L -cosh ˛n x L
		˛1 = 4.694
		˛2 = 7.855
		˛n = n + 1 2	(n > 2)
	Clamped-clamped	˛0 = 4.730		sin ˛n x L -sinh ˛n x L + sin ˛n -sinh ˛n cosh ˛n -cos ˛n	cos ˛n x L -cosh ˛n x L
		˛1 = 7.853
		˛2 = 11.00

[START_REF]OOFELIE::Multiphysics©, User Manual[END_REF] Here again, we try to capture the length of the prongs in the source strength, through what Sillitto calls an « effective length ».