Tao Li

Claude-Henri Lamarque
email: claude.lamarque@entpe.fr

Sébastien Seguy
email: sebastien.seguy@insa-toulouse.fr

Alain Berlioz
email: alain.berlioz@univ-tlse3.fr

Chaotic characteristic of a linear oscillator coupled with vibro-impact nonlinear energy sink

Keywords: Chaos, Lyapunov exponent, Vibro impact, Nonlinear energy sink, Targeted energy transfer, Impact damper

The chaotic characteristic of a system with vibroimpact nonlinear energy sink is studied here. An analytical method is developed to calculate Lyapunov exponent. The mechanism by which impact results in chaos is further clarified rather only the calculation of Lyapunov exponent. In addition, an approach to identify Lyapunov exponents from experimental data is proposed, and the estimated results are consistent with numerical results.

Introduction

Nonlinear systems with impact have been extensively studied for decades [START_REF] Ibrahim | Vibro-impact dynamics: modeling, mapping and applications[END_REF]. Recently, vibro-impact phenomenon has been re-examined from the viewpoint of targeted energy transfer [START_REF] Vakakis | Nonlinear targeted energy transfer in mechanical and structural systems[END_REF]. As a result, its dynamics is further investigated with the application of new analytical methods [START_REF] Lee | Periodic orbits, damped transitions and targeted energy transfers in oscillators with vibro-impact attachments[END_REF][START_REF] Gendelman | Analytic treatment of a system with a vibroimpact nonlinear energy sink[END_REF][START_REF] Gendelman | Dynamics of forced system with vibro-impact energy sink[END_REF][START_REF] Gendelman | Forced system with vibro-impact energy sink: Chaotic strongly modulated responses[END_REF] and even experimental observations [START_REF] Gourc | Targeted energy transfer under harmonic forcing with a vibro-impact nonlinear energy sink: Analytical and experimental developments[END_REF][START_REF] Li | Dynamics of cubic and vibro-impact nonlinear energy sink: Analytical, numerical, and experimental analysis[END_REF][START_REF] Li | On the dynamics around targeted energy transfer for vibro-impact nonlinear energy sink[END_REF][START_REF] Li | Optimization mechanism of targeted energy transfer with vibro-impact energy sink under periodic and transient excitation[END_REF][START_REF]Dynamics of two vibro-impact nonlinear energy sinks in parallel under periodic and transient excitations[END_REF][START_REF] Pennisi | Experimental investigation and analytical description of a vibro-impact nes coupled to a single-degree-of-freedom linear oscillator harmonically forced[END_REF]. Impact is proved to be efficient to control vibration and its corresponding device is termed as vibro-impact (VI) nonlinear energy sink (NES), meanwhile the existence of impact can result in complicated dynamics. For example, the strongly modulated characteristic is experimentally observed from the response of acceleration in [START_REF] Li | On the dynamics around targeted energy transfer for vibro-impact nonlinear energy sink[END_REF] and also observed from the response of displacement in [START_REF] Li | Optimization mechanism of targeted energy transfer with vibro-impact energy sink under periodic and transient excitation[END_REF]. The intermittent feature is well demonstrated by the experimentally observed impact time difference [START_REF] Li | On the dynamics around targeted energy transfer for vibro-impact nonlinear energy sink[END_REF]. Although many aspects of dynamics are further studied, whether the above mentioned characteristic is related to chaos is not verified.

In [START_REF] Gendelman | Dynamics of forced system with vibro-impact energy sink[END_REF][START_REF] Gendelman | Forced system with vibro-impact energy sink: Chaotic strongly modulated responses[END_REF], the strongly modulated response (SMR) of a linear system coupled with VI NES is proved to be chaotic. Actually, the characteristics of modulation and chaos are two different aspects of a single response. The mentioned study is an attempt to combine these two aspects. From the viewpoint of comparing different characteristics, this attempt is pioneering, but this study is not complete for two reasons. Firstly, only the results of Lyapunov exponents are showed and the method to calculate them is not displayed. Secondly, Lyapunov exponent is only an average measure of the divergent or convergent characteristic of two initial close orbits. The relation between impact and chaos should be further explained with more details. Since it is the addition of impact that makes such features possible, the common role of impact that results in different features should be found.

Impact phenomenon is a typical member of nonlinear phenomenon with discontinuity. In the aspect of Lyapunov exponent, a general method to calculate it for dynamic systems with discontinuities is developed by Muller [START_REF] Müller | Calculation of lyapunov exponents for dynamic systems with discontinuities[END_REF] and the transfer matrix at the instants of impact is corrected. With a similar idea as Muller, the chaos of a forced pendulum with friction is studied in [START_REF] Lamarque | Numerical study of a forced pendulum with friction[END_REF]. Around chaos resulting from impact, there are already many studies [START_REF] Stefanski | Estimation of the largest lyapunov exponent in systems with impacts[END_REF][START_REF] Stefański | Estimation of the dominant lyapunov exponent of non-smooth systems on the basis of maps synchronization[END_REF][START_REF] De Souza | Calculation of lyapunov exponents in systems with impacts[END_REF][START_REF] Jin | A method for calculating the spectrum of lyapunov exponents by local maps in non-smooth impactvibrating systems[END_REF]. Their common objective is to calculate Lyapunov exponent in order to judge a response is chaos or not, as has been done for other nonlinear systems, though their methods are different. However, the objective of this paper is not only to calculate Lyapunov exponent, but also to find the factor that contributes to positive Lyapunov exponent. Evidently, it comes from impact, but the specific mechanism is not evident and should be further explored. More specifically, it is to find a more detailed measure of chaos rather than only such a common measure as Lyapunov exponent. That is to say, the role of impact should be directly observed, if chaos is caused by impact. Similarly, the role of friction should be observed if friction is the cause of chaos. In addition, the possibility of identification of Lyapunov exponents from experimental data will be explored here.

The paper is organized as follows. In Section 2, the model of a vibro-impact system is presented. In Section 3, the analytical derivation to calculate Lyapunov exponents is treated. In Section 4, some numerical results are presented. In Section 5, a method is proposed to identify Lyapunov exponents from experimental data. Finally, conclusion of this paper summarizes the main results.

Modeling the vibro-impact system

The considered vibro-impact system is displayed in Fig. 1 [START_REF] Gourc | Targeted energy transfer under harmonic forcing with a vibro-impact nonlinear energy sink: Analytical and experimental developments[END_REF][START_REF] Li | On the dynamics around targeted energy transfer for vibro-impact nonlinear energy sink[END_REF][START_REF] Li | Optimization mechanism of targeted energy transfer with vibro-impact energy sink under periodic and transient excitation[END_REF]. A linear oscillator (LO) is periodically driven by a shaker. A VI NES is coupled with LO only through impact, and the friction between them is neglected.

During periods without impacts, the system is governed by the following equation:

ẍ + ελ ẋ + x = εG sin Ω τ + ε 2 λ GΩ cos Ω τ ε ÿ = 0 ∀ |x -y| < b (1)
and corresponding parameters are as follows:

ε = m 2 m 1 , ω 0 2 = k 1 m 1 , f 0 = ω 0 2π , τ = ω 0 t, λ = c 1 m 2 ω 0 , Ω = ω ω 0 , G = F ε
where x, m 1 , c 1 and k 1 are the displacement, mass, damping and stiffness of the LO, respectively. y and m 2 are displacement and mass of VI NES, respectively. The dots denote differentiation with respect to dimensionless time τ. b represents the clearance. x e is the displacement imposed on the base by shaker. εG sin Ω τ and ε 2 λ 1 GΩ cos Ω τ represent the contribution of force by displacement and velocity of shaker, respectively.

Eq. ( 1) can also be written as:

ẍ + ελ ẋ + x = ε A sin (Ω τ + φ ) ε ÿ = 0 ∀ |x -y| < b (2)
where

A = G 2 + λ 2 ε 2 G 2 Ω 2 φ = arctan (λ ε Ω ) (3)
When |x -y| = b, impacts occur. The relation between after and before impact is obtained under the hypothesis of the simplified shock theory and the condition of total momentum conservation:

x + = x -, y + = y - ẋ+ + ε ẏ+ = ẋ-+ ε ẏ-, ẋ+ -ẏ+ = -R ẋ--ẏ-, for |x -y| = b (4)
where R is the restitution coefficient, and the superscripts + andindicate time immediately after and before impact.

Analytical treatment to calculate Lyapunov exponent

The objective here is to derive formula to calculate the Lyapunov exponents. Correspondingly, it is to calculate the transfer matrix between impacts and the transfer matrix at the moments of impacts.

The state variables of flow Φ can be expressed as:

Φ (τ 0 , x 0 , τ) =     x 1 (τ) x 2 (τ) x 3 (τ) x 4 (τ)     =     ẋ x ẏ y     (5)
then,

Φ = f (Φ) =     -ελ -1 0 0 1 0 0 0 0 0 0 0 0 0 1 0     Φ +     ε A sin (Ω τ + φ ) 0 0 0     (6)
If initial conditions are given:

Φ (τ 0 , x 0 , τ = 0) =     x 1 (0) x 2 (0) x 3 (0) x 4 (0)     =     x10 x20 x30 x40     (7)
The equation of motion between impacts can be obtained and expressed as follows:

x 1 (τ) = ẋ2 (τ) x 2 (τ) = e -ε λ (τ-τ 0 ) (p1 sin (ω d (τ -τ 0 )) + q1 cos (ω d (τ -τ 0 ))) + f1 sin (Ω (τ -τ 0 )) + f2 cos (Ω (τ -τ 0 )) x 3 (τ) = ẋ4 (τ) x 4 (τ) = x30(τ -τ 0 ) + x40 (8)
where τ 0 is the initial time and

f1 = - ε A Ω 2 -1 cos (φ ) -2 λ ε Ω sin (φ ) Ω 4 + (4 ε 2 λ 2 -2) Ω 2 + 1 f2 = - ε A Ω 2 -1 sin (φ ) + 2 λ ε cos (φ ) Ω Ω 4 + (4 ε 2 λ 2 -2) Ω 2 + 1 (9)
and

q1 = x20 -f2 p1 = - -x10 -ε λ x20 + ε λ f2 + f1 Ω ω d (10)
and

ω d = 1 -ε 2 λ 2 (11)
Now, the i th impact is studied, i = 1, 2, 3 • • • . The corresponding time is τ i with τ - i and τ + i denoting the moment before and after this impact.

(τ i ) ≤ τ < τ i+1 Φ = f (Φ), Φ (τ i ) = Φ i (12)
To calculate Lyapunov exponents, an initial difference of states is needed:

𝜏 𝑖+1 𝜏 𝑖+1 = 𝜏 𝑖+1 +𝛿𝜏 𝑥 2 𝜏 -𝑥 4 𝜏 = 𝑏 𝑥 2 𝜏 -𝑥 4 𝜏 = 𝑏 𝛿Φ - 𝛿Φ + Φ (𝜏) Φ(𝜏)
Fig. 2 Trajectory difference before impacts (δ Φ -) and after impacts (δ Φ + ) : original impact in black rectangle (thin broken line) at time τ i+1 in red rectangle (thick broken line), perturbed impact in black rectangle (thin line) at time τ i+1 in red rectangle (thick broken line)

δ Φ(τ) = Φ(τ) -Φ(τ) =     δ x 1 (τ) δ x 2 (τ) δ x 3 (τ) δ x 4 (τ)     (13)
where

Φ(τ) =     x 1 (τ) x 2 (τ) x 3 (τ) x 4 (τ)     , Φ(τ) =     x 1 (τ) x 2 (τ) x 3 (τ) x 4 (τ)     (14)
are state variables of perturbed trajectory and original trajectory, respectively. For the impact of the original motion (τ i+1 ) and that of perturbed motion ( τ i+1 ), the whole process is demonstrated in Fig. 2. The objective here is to calculate the transfer matrix between the trajectory difference at the moment just before this first impact (τ i+1 ) and that just after the moment of the second impact ( τ i+1 ). δ Φ(τ + i ) after the i th impact is:

δ Φ(τ + i ) =     δ x 1 τ + i δ x 2 τ + i δ x 3 τ + i δ x 4 τ + i     (15)
For the period between the i th and (i + 1) th impact, the evolution of differences is governed and can be obtained by the following equation:

(τ i ) ≤ τ < τ i+1 : δ Φ = F (τ) δ Φ (16)
where the starting time is τ + i . With Eq. ( 6),

F (τ) = ∂ f ∂ Φ T | Φ=Φ(τ) =     -ελ -1 0 0 1 0 0 0 0 0 0 0 0 0 1 0     (17)
With Eqs. ( 16) and ( 17),

δ Φ (τ) = J1 i (τ)δ Φ τ + i ( 18
)
where

J1 i (τ) =     a 11 a 12 0 0 a 21 a 22 0 0 0 0 1 0 0 0 τ -τ i 1     (19)
and

a 11 = e -ε λ (τ-τ i ) [cos (ω d (τ -τ i )) - ε λ ω d sin (ω d (τ -τ i ))] a 12 = - e -ε λ (τ-τ i ) ε 2 λ 2 + ω d 2 ω d sin (ω d (τ -τ i )) a 21 = e -ε λ (τ-τ i ) ω d sin (ω d (τ -τ i )) a 22 = e -ε λ (τ-τ i ) ε λ ω d sin (ω d (τ -τ i )) + cos (ω d (τ -τ i )) (20)
The above four coefficients of matrix are not influenced by excitation.

Therefore, the transfer matrix from the moment (τ + i ) to the moment (τ - i+1 ) can be obtained:

J i = J1 i τ - i+1 (21)
Then the transfer matrix at the moment of impact as showed in Fig. 2 is considered. The impact condition is the same for both original and perturbed motions, and can be expressed as follows:

g(τ) = |x 2 (τ) -x 4 (τ)| -b (22)
Then to continue the calculation, suppose x 2 (τ) > x 4 (τ), it just specifies one approaching direction of impact,

G(Φ) = ∂ g(Φ) ∂ Φ T | Φ=Φ(τ) = 0 1 0 -1 (23)
At the moment of impact τ i+1 , the state variables between τ + i+1 and τ - i+1 meet:

    x 1 τ + i+1 x 2 τ + i+1 x 3 τ + i+1 x 4 τ + i+1     = S     x 1 τ - i+1 x 2 τ - i+1 x 3 τ - i+1 x 4 τ - i+1     (24)
where

S =     1-Rε ε+1 0 ε+Rε ε+1 0 0 1 0 0 1+R ε+1 0 -R+ε ε+1 0 0 0 0 1     (25)
S is supposed to be time invariant and is the same at both sides of impacts. According to the work of Muller [START_REF] Müller | Calculation of lyapunov exponents for dynamic systems with discontinuities[END_REF], the following equation can be obtained:

δ Φ τ + i+1 = M i × δ Φ τ - i+1 (26)
where

M i = S - [S Φ(τ - i+1 ) -Φ(τ + i+1 )]G(τ - i+1 ) G(τ - i+1 ) Φ(τ - i+1 ) (27)
where

G(τ - i+1 ) Φ(τ - i+1 ) = x1(τ - i+1 ) -x3(τ - i+1 ) (28)
If x 2 (τ) ≤ x 4 (τ), Eq. ( 27) still holds. Let

δ Φ τ + i+1 = T i × δ Φ τ + i ( 29
)
where

T i = M i × J i ( 30
) when i = 1, δ Φ τ + i = δ Φ (0).
Here the starting time of calculation is fixed to the first impact moment and it will not influence the ultimate goal of the calculation of Lyapunov exponents.

For a period of time Te, we suppose there are n impacts. At the end of n th impact:

δ Φ τ + n = M n × J n • • • M i × J i • • • M 1 × J 1 × δ Φ (0) (31)
Let M(τ i ) be the transfer matrix of the former i impact:

M(τ i ) = M i × J i • • • M 1 × J 1 (32)
Let u i (i = 1..4) be the eigenvalue of transfer matrix M(τ i ) and λ k (k = 1..4) be the Lyapunov exponent.

λ k = 1 τ i -τ 0 ln(|u i |) (33)
where τ 0 is the starting time. Through the calculation of λ k , chaotic performance can be judged from an average viewpoint.

x(0) 0 ẋ(0) 0 y(0) 95% b ẏ(0) 0 4 Numerical results
Here, numerical results related to the calculation of transient Lyapunov exponent are demonstrated.

According to former studies [START_REF] Li | On the dynamics around targeted energy transfer for vibro-impact nonlinear energy sink[END_REF], response regimes can be classified according to the average number of impacts per cycle, and it is noted as z. Here, six categories are showed and described as follows and almost all of them have been experimentally observed.

Type I represents sparse impact without periods of two impacts per cycle and the distribution of impacts is irregular, namely there is only one impact during many cycles. Type II represents SMR with intermittent periods of two impacts per cycle. There are desynchronized parts and synchronized parts. When z is closer to two, the duration of synchronized parts becomes longer. Type III represents regimes with two symmetrical impacts per cycle. Type IV represents regimes with two asymmetrical impacts per cycle. Type V represents other periodic regimes with integer times of impacts per cycle (e.g., three impacts per cycle), symmetrical or asymmetrical. Type VI represents other response regimes with irregular distribution of impacts. Among them, chaotic behavior have been observed for Types I, II and VI.

To start our calculation, the parameters in Table 1 are used and only b is varied to obtain different response regimes.

Type II: chaotic SMR

Firstly, the results of a SMR with b = 0.022 m is showed. The characteristic of a response can be reflected from different viewpoints. All results here are selected only around the calculation of Lyapunov Exponent, and more information can be found in [START_REF] Li | On the dynamics around targeted energy transfer for vibro-impact nonlinear energy sink[END_REF][START_REF] Li | Optimization mechanism of targeted energy transfer with vibro-impact energy sink under periodic and transient excitation[END_REF].

In Fig. 3, the strongly modulated characteristic is well demonstrated by the envelop of displacement. For every peak of x, LO and VI NES move in desynchronization in the increasing part and in synchronization in the decreasing part, i.e., two impacts per cycle. The impact number per cycle is better displayed by relative displacement in Fig. 4 and is equal to the number of extreme values per cycle. In the sparse parts, the amplitude of x increases with z < 2. In the dense parts, the amplitude of x decreases with z = 2.

To calculate Lyapunov exponents, we just calculate one period of limited time as showed by the two broken red line in Fig. 3, considering that too large value of transfer matrix M in Eq. (32) will exceed the capability of computer. In this period, the transient response has disappeared and it should be enough for the calculation of Lyapunov exponents. The evolution of four Lyapunov exponents is showed in Fig. 5. These four transient Lyapunov exponents are almost constant and are positive, which denotes the chaotic characteristic.

However, the evolution of Lyapunov exponents is not smooth, and sometimes there are transient spikes. They result from the calculation of u i in (33). At these non-smooth points, u i is equal or close to zero. Relevant impact points are probably related to the grazing bifurcation but cannot be justified with evidence here. As a result, the calculation error is greatly increased. Nevertheless, the general trend is right and can still predict the evolution of Lyapunov exponents.

To test the influence of the calculation time, different starting and ending times are chosen as showed in Table 2. Although λ 1 and λ 3 vary a little, their difference are small and their sign are always positive. The small difference seems reasonable since every selected time history of displacement is different.

Then, the eigenvalues of matrix J i in Eq. ( 21) without influence of impacts is calculated, and their absolute values are obtained and showed in Fig. 6. The initial difference of LO will decrease with two values less than one as showed in Fig. 6(a-b). This is quite natural since LO is a dissipative system governed by Eq. ( 6) during these periods without impacts. And it is also natural to observe two unit values for VI NES in Fig. 6(c-d) since it is a conservative system itself.

Then, the eigenvalues of matrix M i in Eq. ( 27) at the impact moments are calculated. From the results showed in Fig. 7, the initial difference of displacement will not be changed during the impact process by themselves. On the contrary, the initial difference of velocity will be decreased by impact with a reduction coefficient R by themselves. Since only the values related to eigenvalues are showed here, the interaction between these four states once impacts are not clear. The latter probably could result in any possible chaos. To further observe the influence of impact, the eigenvalues of matrix T i in Eq. ( 30) is calculated. In this way, the interaction of these four states can be observed. As showed in Fig. 8, these two absolute values of eigenvalues of matrix sometimes become very large and greater than one.

Then we compare these eigenvalues with corresponding time history of displacement as showed in Fig. 9, and only the first eigenvalue is demonstrated in Fig. 9(b). It is seen that the eigenvalue is greater than one in area1 and less than one in area2. Here, area1 and area2 are two typical examples of desynchronization and synchronization. For the synchronization parts, the eigenvalue has almost the same value, but this value changes a lot in the desynchronization parts. In area3, the largest value is observed and it has a close relation to the maximal amplitude of modulated displacement.

Evidently, the chaotic characteristic results from impact, but only the existence of impact will not definitively result in chaos. From the above results, chaos results from irregular impacts, namely during desynchronization between LO and VI NES. Since a Lyapunov exponent reflects an average trend, the results here demonstrate more details and reveal the essence of chaos. About this point, it will be further verified from the following results.

Other types

Type I: the evolution of Lyapunov exponents with b = 0.04 m is showed in Fig. 10, four positive values are obtained. Therefore, the irregular and occasional impacts can result in chaos. 11 and Fig. 12, respectively. The former is with two symmetrical impacts and the latter asymmetrical. For these two, four negative values (two of them almost zero) are obtained. Therefore, they are not chaos and the regular impact does not generate chaotic behavior. Then b is decreased to 0.004 m, three impacts per cycle is obtained, namely Type V. The evolution of eigenvalues is showed in Fig. 13 and four negative values are also obtained. When b is further increased, the impact number is between 3 and 4, as showed in Fig. 14(a). For this irregular response, the final Lyapunov exponents are positive. 3. With the increase of b, the duration of two impacts per cycle decreases and transient Lyapunov exponents increase, which verifies the initial assumption.

Experimental identification of Lyapunov Exponents

In [START_REF] Li | On the dynamics around targeted energy transfer for vibro-impact nonlinear energy sink[END_REF][START_REF] Li | Optimization mechanism of targeted energy transfer with vibro-impact energy sink under periodic and transient excitation[END_REF], the above-mentioned response regimes have been experimentally observed. Here, the objective is to explore the possibility of the identification of Lyapunov exponents from experimental data. Firstly, a theoretical method to identify Lyapunov exponents from experimental data is presented. Then, it is applied to both a steady response with two impacts per cycle and a SMR with chaotic characteristics.

Analytical development

Based on the analytical development in Section 3, we propose a method here to identify Lyapunov exponents from experimental data in Ref. [START_REF] Li | On the dynamics around targeted energy transfer for vibro-impact nonlinear energy sink[END_REF].

To calculate Lyapunov exponents of a specific response, matrix J i and M i should be estimated from experimental data. The common parameters in Table 4 can be used as the starting point of calculation, and other parameters should be identified from experimental data and they are further calculated to finish the estimation process. Specifically, to calculate J i , only the impact time are further needed. Then to calculate M i in Eq. ( 27), the situation will be a little complex and some approximations will be used. Among the parameters needed to calculate M i , G is constant as follows:

G = 0 1 0 -1 (34)
S(τ) will be a time variable if R is not time invariant or has different values at different impact sides.

S(τ) =     1-Rε ε+1 0 ε+Rε ε+1 0 0 1 0 0 1+R ε+1 0 -R+ε ε+1 0 0 0 0 1     (35)
M i is expressed in the following form:

M i = S(τ i+1 ) - [S(τ i+1 ) Φ(τ - i+1 ) -Φ(τ + i+1 )]G G Φ(τ - i+1 ) (36)
Because ÿ(τ) = 0 , the state variables before and after impact are:

Φ(τ - i+1 ) = ẍ(τ - i+1 ) ẋ(τ - i+1 ) 0 ẏ(τ - i+1 ) (37)
Φ(τ + i+1 ) = ẍ(τ + i+1 ) ẋ(τ + i+1 ) 0 ẏ(τ + i+1 ) (38)
and

G Φ(τ - i+1 ) = ẋ(τ - i+1 ) -ẏ(τ - i+1 ) (39)
If x, ẋ and ẍ are given, y should meet the following requirement at any impact moments:

|x(τ i ) -y(τ i )| = b (40)
If the direction of relative displacement between LO and VI NES is changed, the following relation should hold:

(x(τ i ) -y(τ i ))(x(τ i+1 ) -y(τ i+1 )) < 0 (41)
The reason why Eq. ( 40) is simplified to the form in Eq. ( 41) is that almost all relative displacement change direction in our experimental observations. If other cases are encountered, Eq. (40) can still be applied.

Supposing that the friction between VI NES and LO is small enough that it will not change the velocity of VI NES, then the following approximation relation between y and ẏi can be obtained:

y(τ i+1 ) -y(τ i ) = ẏi (τ i+1 -τ i ) ( 42
)
where ẏi is the velocity at i th time interval [τ i , τ i+1 ], and combining Eqs. (40) and (42), ẏi can be calculated:

ẏi = x(τ i+1 ) -x(τ i ) ± 2b τ i+1 -τ i (43)
With obtained ẏ, restitution coefficient R can be calculated by Eq. (4).

Application

To estimate Lyapunov exponents, an estimated constant restitution coefficient with R = 0.85 is applied here, and this value of restitution coefficient is just an approximation to its real value and is just enough to demonstrate the proposed method here.

Non-chaotic

For b=9mm, a steady response with two impacts per cycle is observed and its time history of acceleration is showed in Fig. 15, in which two impacts per cycle is well demonstrated.

Fig. 15 Acceleration of LO with b = 0.009 m: impact moments are denoted by cross [START_REF] Li | Dynamics of cubic and vibro-impact nonlinear energy sink: Analytical, numerical, and experimental analysis[END_REF] Its time history of displacement is showed in Fig. 16 and impacts moments are marked out by red circles. Then the velocities before and after impacts are estimated, and the corresponding phase portrait between displacement and velocity is showed in Fig. 17 Because these values of displacements and velocities after and before impacts are so close that they cannot be clearly distinguished. To resolve this problem, the displacements and velocities after impact at the Poincare section with |x -y| = b is showed in Fig. 18.

Then Lyapunov exponents are estimated during limited time with limited impact numbers, and the result is showed 11 and 12, it can be anticipated that it will approach zero so long as the duration of time is long enough.

Chaotic

For b=15mm, a chaotic response is observed. Its chaotic characteristic can be observed from different viewpoints. Its time history of acceleration is showed in Fig. 20, and a detailed analysis of its chaotic impact time difference can be found in [START_REF] Li | Dynamics of cubic and vibro-impact nonlinear energy sink: Analytical, numerical, and experimental analysis[END_REF].

Its time history of displacement is showed in Fig. 21 and impacts moments are marked out by red circles.

Then, the velocity is estimated from displacement and acceleration, and the corresponding phase portrait between displacement and velocity is showed in Fig. 22. The chaotic Moreover, it is better demonstrated by the displacements and velocities after impacts at the Poincare section with |x -y| = b as showed in Fig. 23. Scattered points at both places are a well proof of chaos compared to these of the former periodic response in Fig. 18.

Then Lyapunov exponents are estimated during this limited time with limited impact numbers, and the result is showed in Fig. 24. Two facts are observed as follows: the transient Lyapunov exponents are positive and they are almost already stable. Compared to the results in Fig. 5, they are consistent.

It has to be pointed out that the same conclusion can be obtained for other restitution coefficients for the impacts at the two sides of LO or even different restitution coefficients at these two sides. In this paper, the method to calculate Lypapunov exponents of a vibro-impact system is theoretically derived. Combined with numerical simulation, the chaotic characteristics of different responses are studied. Finally, an approach to identify Lyapunov exponents from experimental data is proposed and is analyzed with some existing experimental data.

Lyapunov exponents for different response regimes is obtained and the chaotic characteristic is verified. In addition, it is found that the essence of chaotic is the irregular and desynchronized impacts, which is common base of the same chaotic feature of different response regimes. Therefore, although the overall effect reflected by Lyapunov exponents is important, the underlying factor that contributes to this average effect should be found. By this way, the relation between chaotic characteristic and other characteristic such as modulation and synchronization could be clearer. This point is more fundamental than Lyapunov exponents. This idea applies for other nonlinear systems with chaos.

Then, Lyapunov exponents are estimated for both a steady response with two impact per cycle and SMR, i.e. non-chaotic and chaotic response. Transient positive Lyapunov exponents within a limited time are obtained for both responses. A continuously decreasing positive value is obtained for nonchaotic response in this limited time and a steady positive value is obtained for SMR. Therefore, this response regime can be termed as chaotic SMR. The experimental results are consistent with what have been numerically obtained. But it has to be pointed out that the restitution coefficient will not be constant from an experimental viewpoint, and the used fixed value here is chosen as close to the fact as possible. The main error comes from different impact materials at two sides of LO and the neglected friction between VI NES and LO. The specific way to estimate restitution coefficient is out the scope of this paper and will be presented later.

Fig. 1

 1 Fig. 1 Schema of a LO coupled with a VI NES under periodic excitation

Fig. 3

 3 Fig. 3 Type II and SMR with b = 0.022 m: displacement x of LO

Fig. 5 Fig. 7

 57 Fig. 5 Type II and SMR with b = 0.022 m: evolution of four Lyapunov exponents

Fig. 8

 8 Fig. 8 Type II and SMR with b = 0.022 m: absolute value of four eigenvalues of matrix T i

Fig. 10

 10 Fig. 10 Type I with b = 0.04 m: evolution of Lyapunov exponents. (a) Relative displacement yx; (b) transient Lyapunov exponents

Fig. 11

 11 Fig. 11 Type III with b = 0.015 m: evolution of Lyapunov exponents. (a) Relative displacement yx; (b) transient Lyapunov exponents

Fig. 13

 13 Fig. 13 Type V with b = 0.004 m: evolution of Lyapunov exponents. (a) Relative displacement yx; (b) transient Lyapunov exponents

Fig. 14

 14 Fig. 14 Type VI with b = 0.003 m: evolution of Lyapunov exponents. (a) Relative displacement yx; (b) transient Lyapunov exponents

  .

Fig. 16

 16 Fig. 16 Displacement of LO with b = 0.009 m: impact moments are denoted by circles

Fig. 17

 17 Fig. 17Phase portrait of LO between displacement and velocity with b = 0.09 m: the moments before and after impact are denoted by red circles and crosses

  Fig. 17Phase portrait of LO between displacement and velocity with b = 0.09 m: the moments before and after impact are denoted by red circles and crosses

Fig. 19

 19 Fig. 19 Evolution of Lyapunov exponents of a steady response with two impacts per cycle for b = 0.009 m

Fig. 20

 20 Fig. 20 Acceleration of LO with b = 0.015 m: impact moments are denoted by cross[START_REF] Li | Dynamics of cubic and vibro-impact nonlinear energy sink: Analytical, numerical, and experimental analysis[END_REF]

Fig. 21

 21 Fig. 21 Displacement of LO with b = 0.015 m: impact moments are denoted by circles

Fig. 22

 22 Fig. 22 Phase portrait of LO between displacement and velocity with b = 0.015 m: the moments before and after impact are denoted by red circles and crosses

Fig. 23

 23 Fig. 23 Poincare section with |x -b| = b: two groups of velocities and displacement of LO after impact

Fig. 24

 24 Fig. 24 Evolution of Lyapunov exponents of a SMR with b = 0.015 m

Table 1

 1 Simulation parameters

		Physical Parameters
	m 1	3.807 kg	c 1	2.53 Ns/m
	k 1	11.68 * 10 3 N/m	m 2	32 g
	R	0.6
		Reduced Parameters
	ε	0.84 %	λ	1.43
		Excitation Parameters
	G	0.02	Ω	1
		Initial conditions

Table 2

 2 Transient Lyapunov exponents at different time periods [τ s , τ e ]

		λ 1	λ 3
	[18.1300 1708.9] 0.1093 0.08427
	[1708.9 3652.7]	0.1143 0.09242
	[3652.7 5521.3]	0.1197	0.1015

Table 3

 3 Transient Lyapunov exponents of different regimes

	b	λ 1	λ 3
		Type I
	0.04	0.1321	0.1245
	0.034	0.1372	0.127
		Type II: SMR
	0.031	0.1528	0.1391
	0.028	0.1452	0.1309
	0.025	0.1138	0.09727
	0.022	0.1093	0.08427
	0.020	0.1014	0.0783
	0.018	0.06368 0.03913
	0.017	0.04703 0.02043
	0.0165 0.03686 0.01084

Table 4

 4 Experimental Parameters[START_REF] Li | On the dynamics around targeted energy transfer for vibro-impact nonlinear energy sink[END_REF]

		Physical Parameters
	m 1	4.168 kg	c 1	3.02 Ns/m
	k 1	11.47 * 10 3 N/m m 2	32 g
		Reduced Parameters
	ε	0.76%	λ	1.80
	of two impacts per cycle should possess different Lyapunov
	exponents. This duration is shorter, positive Lyapunov expo-
	nents should be larger. For this reason, transient Lyapunov
	exponents of different response regimes with different b are
	put in Table

Acknowledgements The authors acknowledge the French Ministry of Science and the Chinese Scholarship Council under grant No. 20130449 0063 for their financial support.


Fábio M Costa

Nikolaos Georgantas
email: nikolaos.georgantas@inria.fr

Raphael De A Gomes
email: raphael.gomes@ifg.edu.br

Ricardo C A Da Rocha
email: rcarocha@acm.org

Georgios Bouloukakis

Cross-layer QoS-Aware Resource Allocation for IoT-Enabled Service Choreographies

Keywords: • Applied computing → Service-oriented architectures, • Software and its engineering → Middleware, Middleware, Service choreographies, QoS, Resource allocation

Resource allocation for distributed systems is a well-known approach to deal with quality of service requirements. However, existing approaches do not consider the effects of resource allocation at the different levels of a system, especially when considering the end-to-end behavior of distributed compositions such as service choreographies. Based on our previous research on the optimized and QoS-aware resource allocation for service choreographies and publish-subscribe brokers, we outline a novel approach that integrates the treatment of resource allocation at both levels in order to yield a more precise control over end-to-end QoS. The paper presents an early architecture and discusses some of the main challenges towards realizing the approach.

INTRODUCTION

Large-scale distributed applications in the Future Internet will involve a multitude of heterogeneous services and resources that will Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s). M4IoT'18, December 10-11, 2018, Rennes, France © 2018 Copyright held by the owner/author(s). ACM ISBN 978-1-4503-6118-7/18/12. https://doi.org/10.1145/3286719.3286726 need to be integrated and coordinated in a dynamic way to match evolving requirements and changing execution environments. Effective coordination in this context can be achieved through service choreography protocols [START_REF] Barker | Choreographing Web Services[END_REF], which enable services to coordinate themselves in a decentralized way to achieve meaningful collaboration. Nevertheless, achieving such distributed coordination in a way that preserves end-to-end QoS and makes efficient use of resources remains an open research problem.

In previous work, we proposed an end-to-end QoS-aware approach for the optimized allocation of resources to the services that take part in choreographies [START_REF] Gomes | A Model-Based Approach for the Pragmatic Deployment of Service Choreographies[END_REF]. While the approach is effective in guaranteeing the QoS of service choreographies, its precision is limited due to the fact that it abstracts away the effect of the middleware layer on end-to-end QoS. Middleware components, such as pub/sub brokers, message queues, and service buses are an integral part of the data path, and thus also contribute to performancerelated QoS attributes, such as end-to-end delay and throughput. A finer level of control over the QoS of service choreographies can thus be achieved by considering the middleware components, together with the application-level services as first-class entities with respect to resource allocation.

While the problem of QoS-constrained resource allocation for service choreographies (and indeed for service compositions in general) is not new, the approach presented in this paper differs from existing approaches by proposing a novel integration of four complementary architectural elements. Firstly, we adopt a cross-layer view of resource allocation, by considering both the applicationand middleware-level entities involved, enabling greater control over the effects of resource allocation on QoS. Secondly, we use a proactive approach, which considers the service demand imposed by choreography clients, together with the capacity and QoS requirements of the services in a end-to-end manner prior to the actual choreography enactment, thus avoiding excessive QoS violations while still allowing runtime adaptation if needed. This is in contrast with reactive approaches, such as in [START_REF] Furtado | Towards an Enactment Engine for Dynamically Reconfigurable and Scalable Choreographies[END_REF] and [START_REF] Bucchiarone | Design for adaptation of servicebased applications: Main issues and requirements[END_REF], which start with a initial (default) resource configuration and adapt it as QoS violations are detected. Thirdly, as introduced in previous work [START_REF] Gomes | QoS-Aware Resource Allocation for Mobile IoT Pub/Sub Systems[END_REF], resource allocation at the middleware level takes into account the possibility of intermittent connectivity of services, thus contributing to make the approach feasible for service choreographies that involve mobile services and IoT devices. Finally, resource allocation takes into account the fact that services can be shared across multiple choreographies, enabling the optimized allocation of resources [START_REF] Gomes | A Model-Based Approach for the Pragmatic Deployment of Service Choreographies[END_REF].

In this position paper, we review our previous work on enabling QoS-aware resource allocation, separately, for the application services layer and for the middleware layer, which in turn provides the building blocks for the integrated cross-layer resource allocation approach that we propose. We discuss a motivating scenario and its handling using the cross-layer approach in Section 2. In Section 3 we describe the foundation elements of the proposed approach, while Section 4 outlines an architectural model to realize it. In Section 5 we present a research agenda in the form of open research questions that need to be answered in order to achieve a concrete realization of the approach. Finally, Section 6 concludes the paper.

MOTIVATION

Consider the scenario depicted in Figure 1, which shows an application for controlling urban traffic based on the processing of video streams from cameras distributed across a city area. This application may be implemented as a choreography of IoT devices and services, without a central point of coordination and typically exploiting the flexible cloud-based allocation of resources for services.

Figure 1: Motivation example of a choreography of services and IoT devices

The diversity and scale in this scenario suggest the adoption of appropriate communication middleware (e.g., pub/sub), in order to deal with heterogeneity, mobility, intermittent connectivity, faults, and scalability. Hence, the communication middleware and its corresponding allocated resources not only impact the effectiveness of message delivery but also its end-to-end QoS, particularly in terms of message delays and throughput.

The allocation of suitable resources for the participating services is a common strategy to improve QoS, and it is relatively wellsupported by cloud providers. However, each service is typically considered in isolation, both from other services and from the underlying middleware, thus lacking coordination when allocating resources for an entire choreography of services, and not allowing the system to consider overall goals in terms of end-to-end QoS. Moreover, dealing with the intermittent connectivity of mobile services and IoT devices may require different resource allocation actions at both levels, e.g., to deal with voluntary disconnections at the application services level, and with forced network disconnections at the middleware level. This motivates the need for jointly managing the allocation of resources for both services and their underlying communication middleware (e.g., publish-subscribe brokers), which in turn may require managing the contribution of each service and middleware component towards the end-to-end QoS.

We organize the resource allocation problem in two layers of QoS concerns, depicted in Figure 2: (i) application layer, which comprises each service or IoT device that participates in the choreography; and (ii) middleware layer, which comprises the components of the middleware infrastructure (e.g., brokers) that enable the communication links among services and IoT devices. We argue for a holistic approach for resource allocation that acts at both layers, taking into account the contribution of each application-or middlewarelevel component towards the end-to-end QoS. Hence, an effective strategy for resource allocation should carry out coordinated crosslayer actions, interrelating local choreography adaptations with their effect on the global, end-to-end, QoS.

EARLY EVIDENCE

From the above discussion, it is clear that accurate QoS-aware resource allocation depends on taking into account both the application and middleware layers. In this section, we outline elements of our previous work that consider each of these levels separately, and which can be used as the basis for a cross-layer approach.

Resource Allocation for Services

In previous work [START_REF] Gomes | A Model-Based Approach for the Pragmatic Deployment of Service Choreographies[END_REF], we addressed the problem at the application layer, by proposing a QoS-aware approach for resource estimation and selection when deploying multiple service choreographies with shared services. Given a high-level description of service choreographies and related constraints, the approach enables autonomic resource allocation in a hybrid cloud environment, with multiple cloud providers, whilst decreasing resource utilization costs.

We formalize the resource allocation problem for service choreographies and solved it using a matching strategy based on the Multiple-choice Multi-dimension Knapsack Problem (MMKP) [START_REF] Moser | An algorithm for the multidimensional multiple-choice knapsack problem[END_REF]. We use a graph structure to represent the combined view of all service choreographies, so as to enable the matching between services and available resources. As part of this process, it is necessary to estimate the QoS contribution of each service on all end-to-end paths in the multiple choreographies in which the service takes part. For simplicity, we only considered performance related QoS attributes. In order to estimate the QoS contribution of each candidate resource we rely on queuing theory [START_REF] Edward D Lazowska | Quantitative system performance: computer system analysis using queueing network models[END_REF].

As illustrated in Figure 3, resource selection is carried out by analyzing each end-to-end interaction flow in the service choreography definitions provided as input. For each pair of service and candidate resource, the expected QoS contribution in the flow is estimated using as input the load from the related choreography, as well as the cumulative load from other choreographies that also use the same service. Moreover, to model the processing rates of messages passing through each service, we use an M/M/c queuec is the number of processor cores of each candidate resource and µ pr is the rate for processing messages at each core.

In [START_REF] Gomes | A Model-Based Approach for the Pragmatic Deployment of Service Choreographies[END_REF], we focused only on application-level services and abstracted away the effects of middleware components on resource utilization as well as on the transmission delay when carrying out the performance analysis. By doing so, we were able to focus on the specificities of service sharing among choreographies. However, we recently extended the proposed approach to target resource allocation for the middleware layer, as described next.

Resource Allocation for Pub/Sub Brokers

Publish/Subscribe middleware protocols decouple peers (publishers and subscribers) in both time and space. They increase successful interactions especially when peers are intermittently connected. Additionally, they introduce mechanisms for event routing and transmission. To support distributed applications spanning a wide-area, the pub/sub system has to be implemented as a set of independent, communicating brokers, forming a broker overlay. An important design-time decision is the deployment of each broker to an appropriate machine (resource allocation). While a number of resource allocation solutions exist, they mainly focus on the satisfaction of local QoS requirements (e.g., employing load balancing strategies for individual jobs) lacking support for end-to-end non-functional properties.

In [START_REF] Gomes | QoS-Aware Resource Allocation for Mobile IoT Pub/Sub Systems[END_REF], we provide a QoS-aware solution for allocating resources to pub/sub brokers by taking into account the end-to-end response time between peers, as well as their intermittent connectivity. As shown in Fig. 4, we follow a similar approach as the one described in Section 3.1 by modeling the middleware layer of interacting peers and using a matching strategy. In particular, by relying on [START_REF] Bouloukakis | Timeliness evaluation of intermittent mobile connectivity over pub/sub systems[END_REF][START_REF] Bouloukakis | Simulation-based Queueing Models for Performance Analysis of IoT Applications[END_REF] we model each broker by using an input and an output queue. The input queue is an M/M/c (µ pr is the core's event processing rate). The output queue can be either an M/M/1 or an ON /OF F queue. In both cases the µ t r is used to model the transmission delay of events inside the network. ON /OF F output queues are used to model intermittently connected subscribers, while M/M/1 queues are employed to model event transmission between brokers.

It is worth noting that in [START_REF] Gomes | QoS-Aware Resource Allocation for Mobile IoT Pub/Sub Systems[END_REF] we ignore possible application layer event processing, which is targeted in the approach proposed here.

OVERALL ARCHITECTURE

In our ongoing work, we are combining the building blocks described in the previous section as part of an integrated approach for cross-layer QoS-aware resource allocation for service choreographies. We are currently considering choreographies where services are connected by pub/sub middleware, notably aiming at supporting services that represent IoT devices. We now outline the main elements of the system architecture.

Our initial proposal for the combined approach is to model both levels as a single optimization problem where the goal is to jointly allocate resources for application-level services and for message brokers, subject to common QoS constraints. In doing so, the joint approach matches the system elements (application-level services and middleware components) to resources taking into the intermittent behavior of communicating peers. As illustrated in Figure 5, the input for the proposed matching strategy includes:

• the system model, which comprises a choreography-like description of the service interaction flows and inter-service communication paths, as well as the queuing models described in the previous section and the related QoS constraints; and • the resource model, which describes the available resource types in a multi-cloud environment.

The aforementioned approach relies on the fact that the individual solutions presented in Section 3 are essentially based on the same strategy. More precisely, we abstract the application-and middleware-level entities to be deployed/inspected as a group where each entity can participate in multiple contexts, such as multiple service choreographies in the application layer and multiple message routing paths in the middleware layer. However, some open research questions must be answered (as highlighted next), especially if we take into account heterogeneous communication paradigms instead of only pub/sub.

OPEN RESEARCH QUESTIONS

An initial concrete realization of the approach proposed in this paper part of our ongoing work. We are currently investigating the interaction of the queuing models used at the two levels (application and middleware) considering end-to-end paths across multiple service choreographies and the use of pub/sub middleware. Although in previous work we have validated the queuing network model for arbitrary paths on a network of pub/sub brokers, we relied on predefined paths between the communicating entities. Providing a solution that considers the multiple, possibly dynamic paths (and their corresponding loads) between services in a multi-choreography scenario is still an open issue.

We also need to refine the system model to represent the mapping between choreography flows and communication paths at the middleware level. The model should consider the possibility of having multi-hop paths at the middleware level to connect adjacent services in a choreography. Moreover, it must cater for the possibility that these paths may change dynamically.

Another important issue relates to the likelihood that multiple middleware technologies, representing different interaction paradigms (including client-server middleware and message queues), may be needed in large-scale distributed service choreographies. The QoS and resource allocation models need to consider such heterogeneity. For instance, while the middleware layer in the case of pub/sub brokers is typically allocated on separate hosts, for clientserver, the middleware components are usually collocated with the application components (e.g., in the form of client and server libraries). Moreover, it may also be necessary to consider the resources requirements of the interoperability mechanisms used to connect the different middleware paradigms.

Finally, while the proposed approach assumes that applicationlevel services and middleware components can be seamlessly treated as first class entities, further research is needed to refine this assumption. Notably, while application services can be shared across multiple choreographies, middleware components can be shared by different services that take part on the same or different choreographies. That is, multiplexing is more complex at the middleware level, and we need to investigate how to address this issue in the system model. Similarly, regarding the treatment of intermittent connectivity proposed in [START_REF] Gomes | QoS-Aware Resource Allocation for Mobile IoT Pub/Sub Systems[END_REF], and as suggested in Section 2, we need to investigate whether a simple extension of the use of ON/OFF queues to the application level (to model the voluntary disconnections of mobile services and IoT devices) would be consistent with the treatment of disconnections at the middleware level, especially on an end-to-end basis.

CONCLUSION

In this position paper, we leverage on previous research to propose a novel approach for QoS-aware resource allocation for the application and middleware components that take part in distributed service choreographies. The approach takes into account the possibility of having services with intermittent connectivity, thus paving the way to support choreographies that involve IoT devices and mobile services. The existing building blocks were described, followed by an initial architecture that integrates them. We are currently working on a concrete realization in order to validate the approach and shed light on the research challenges identified so far.

Figure 2 :

 2 Figure 2: Cross-layer resource allocation

Figure 3 :

 3 Figure 3: Performance model for the application layer in service choreographies.

Figure 4 :

 4 Figure 4: Performance model for the middleware layer in pub/sub systems.

Figure 5 :

 5 Figure 5: Overall architecture.

ACKNOWLEDGMENTS

The authors would like to thank FAPEG and Inria for funding the ACHOR Associate Team, within which this work was developed.


H Boujlida
email: boujlida.hanen12@gmail.com

H Haddar †

M Khenissi
email: moez.khenissi@gmail.com

THE ASYMPTOTIC OF TRANSMISSION EIGENVALUES FOR A DOMAIN WITH A THIN COATING

Keywords: transmission eigenvalues, asymptotic expansions, thin layers, inverse scattering problems AMS subject classifications. 35P30, 35P25

We consider the transmission eigenvalue problem for a medium surrounded by a thin layer of inhomogeneous material with different refractive index. We derive explicit asymptotic expansion for the transmission eigenvalues with respect to the thickness of the thin layer. We prove error estimate for the asymptotic expansion up to order 1. This expansion can be used to obtain explicit expressions for constant index of refraction.

1. Introduction. This work is a contribution to the study of transmission eigenvalues [START_REF] Colton | Inverse Acoustic and Eletromagnetic Scattering Theory[END_REF][START_REF] Cakoni | Haddar Inverse Scattering Theory and Transmission Eigenvalues SIAM publications[END_REF][START_REF] Cakoni | Transmission eigenvalues in inverse scattering theory Inverse Problems and Applications[END_REF] and their relation to the shape and material properties of scatterers. These special frequencies are associated with the existence of an incident field that does not scatter. They can be equivalently defined as the eigenvalues of a system of two coupled partial differential equations posed on the inclusion domain. One of these equations refers to the equation satisfied by the total field and the other one is satisfied by the incident field. The two equations are coupled on the boundary by imposing that the Cauchy data coincide. This eigenvalue problem can then be formulated as an eigenvalue problem for a non-selfadjoint compact operator. Although non intuitive, it can be shown that this problem admits an infinite discrete set of real eigenvalues without finite accumulation points [START_REF] Cakoni | The existence of an infinite discrete set of transmission eigenvalues[END_REF][START_REF] Pavarinta | Transmission eigenvalues[END_REF]. These special frequencies can be identified from far field data as proved in [START_REF] Cakoni | On the determination of dirichlet or transmission eigenvalues from farfield data[END_REF][START_REF] Kirsch | The inside-outside duality for scattering problems by inhomogeneous media[END_REF][START_REF] Cakoni | Haddar Inverse Scattering Theory and Transmission Eigenvalues SIAM publications[END_REF]. Since they carry information on the material properties of the scatterer, transmission eigenvalues would then be of interest for the inverse problem of retrieving qualitative information on the material properties from measured multistatic data [START_REF] Giorgi | Computing estimates of material properties from transmission eigenvalues[END_REF][START_REF] Harris | Transmission eigenvalues and non-destructive testing of anisotropic magnetic materials with voids[END_REF]. In this perspective, it appears important to study the dependence of these eigenfrequencies with respect to the material properties and the geometry. Several works in the literature have addressed this issue by considering asymptotic regimes and quantifying the dependence of the first leading terms in the asymptotic expansion of the transmission eigenvalue with respect to the small parameter [START_REF] Cakoni | Asymptotic analysis of transmission eigenvalue problem for a Dirichlet obstacle coated by a thin layer of a non-absorbing media[END_REF][START_REF] Cakoni | Asymptotic expansions for transmission eigenvalues for media with small inhomogeneities[END_REF][START_REF] Moskow | Nonlinear eigenvalue approximation for compacts operators[END_REF][START_REF] Cakoni | The imaging of small perturbations in an anisotropic media[END_REF]. We here consider the case of a scatterer made of a thin coating which corresponds to frequently encountered configurations in the stealth technology for instance. The goal is to characterize the dependence of the first order term on the material properties and the thickness of the coating. A first work on this topic was done in [START_REF] Cakoni | Asymptotic analysis of transmission eigenvalue problem for a Dirichlet obstacle coated by a thin layer of a non-absorbing media[END_REF] where the case of coated perfect scatterer is considered.

One proves in particular for the latter case that the first order term depends only on the thickness.

We here address the more complicated configuration of a coated penetrable media. The analysis indicates that the first order asymptotic resembles to the shape derivative for the buckling plate equation [START_REF] Henrot | Extremum Problems for Eigenvalues of Elliptic Operators[END_REF] and contain non trivial dependence on the material properties. More importantly, this expansion allows us to obtain explicit (approximate) expressions for the thin layer index of refraction in terms of the thickness of the layer, the transmission eigenvalue for the coated medium that can be extracted from the measurements and the transmission eigenvalues and eigenvectors for the coated free medium that can be evaluated numerically. This indeed can be useful for the solution of the inverse problem.

Although the formal derivation follows the systematic procedure using the classical scaled expansion method (as in [START_REF] Bendali | The effect of a thin coating on the scattering of a time-harmonic wave for the Helmholtz equation[END_REF][START_REF] Aslanyürek | Generalized impedance boundary conditions for thin dielectric coatings with variable thickness[END_REF][START_REF] Delourme | Approximate models for wave propagation across thin periodic interfaces[END_REF] for instance), the rigorous justification is much more involved.

For instance the arguments in [START_REF] Cakoni | Asymptotic analysis of transmission eigenvalue problem for a Dirichlet obstacle coated by a thin layer of a non-absorbing media[END_REF] are hard to extend to the present case since special uniform estimates have to be obtained for the transmission problem. We restrict ourselves here to the justification of the first two terms in the asymptotic expansion using the abstract theory developed in [START_REF] Je Osborn | Spectral approximations for compacts operators[END_REF][START_REF] Moskow | Nonlinear eigenvalue approximation for compacts operators[END_REF]. We follow the procedure developed in [START_REF] Cakoni | Asymptotic expansions for transmission eigenvalues for media with small inhomogeneities[END_REF] for the case of small obstacles asymptotic.

The main technical point in the proof is to obtain the corrector for the main operator, which is here the biharmonic operator. Our main result provides explicit expansion for simple transmission eigenvalues and for multiple transmission eigenvalues that are associated with a generalized eigenspace spanned only by eigenvectors.

We analyze the problem where the contrast in material properties affect only the lower order term in the Helmholtz equation. We finally indicate that although the problem is considered only in dimension 2, the results of the main theorem (including the expression of the first order asymptotic term) remain true for three dimensions (up to more complicated technicalities in the proof related to differential geometry).

The paper is organized as follows. We first introduce the transmission eigenvalues and write them as the eigenvalues of a non selfadjoint operator. We then present the main result of our paper and discuss applications to the inverse problem. We present next the outline of a classical formal procedure to obtain the expression of the asymptotic expansion. We give the expression till the second order term. We explain in particular why the expression of the second order term would have less interest in practice. We then proceed with the main part of the paper that provides explicit expressions and an error estimate for the first two terms in the asymptotic expansion.

Problem statement and main results.

Let Ω ⊂ R 2 be a bounded domain with a smooth boundary Γ. We denote by

Ω 0 = {x ∈ Ω, d(x, Γ) > } and its boundary Γ = {x ∈ Ω, d(x, Γ) = } = ∂Ω 0 ,
for > 0 a small enough parameter, where d(x, Γ) denotes the distance of a point x to the boundary Γ. Let Ω = Ω\Ω 0 be the layer of thickness around Ω 0 (see Figure 1). We consider the following transmission eigenvalue problem:

(1)

         ∆w + k 2 n (x)w = 0 in Ω, ∆v + k 2 v = 0 in Ω, w = v on Γ, ∂w ∂ν = ∂v ∂ν on Γ,
where k denotes the unknown eigenfrequency and ν the unitary normal to Γ directed to the interior of Ω. The index of refraction n is defined as follows:

n (x) = n 0 (x) in Ω 0 , n 1 (x) in Ω ,
where n 0 and n 1 are non negative real valued functions ∈ L ∞ (R 2 ) that are independent from .

For the sake of simplicity, we assume that the restriction of n 0 and n 1 to Ω are constant functions along the normal coordinate to Γ for sufficiently small. We finally assume that the function 1/(1 -n ) is either positive definite or negative definite on Ω. We remark that this assumption also implies that 1/(1 -n 0 ) is either positive definite or negative definite on Ω and that

(2) 1/|1 -n (x)| ≥ γ > 0 for a.e.
x ∈ Ω with γ being independent from (sufficiently small) .

The main goal of this paper is to find the asymptotic expansion of eigenfrequencies k with respect to . Assuming that 1 1 -n ∈ L ∞ (Ω), the transmission eigenvalue problem (1) can be reformulated as the nonlinear eigenvalue problem for λ := k 2 ∈ R and u := w -v ∈ H 2 0 (Ω) such that

(∆ + λ n ) 1 1 -n (∆ + λ )u = 0 in Ω,
which in variational form, after integration by parts, is formulated as finding λ ∈ R and non-trivial

function u ∈ H 2 0 (Ω) such that (3) Ω 1 1 -n (∆u + λ u )(∆φ + λ n φ)dx = 0, ∀φ ∈ H 2 0 (Ω).
The space H 2 0 (Ω) denotes the closure in H 2 (Ω) of the set of regular compactly supported functions in Ω. We shall work with the reformulation of (3) as a linear eigenvalue problem for a non selfadjoint compact operator [START_REF] Cakoni | Haddar Inverse Scattering Theory and Transmission Eigenvalues SIAM publications[END_REF]. First observe that (3) can be written as

(4) A u + λ B u + λ 2 C u = 0 in H 2 0 (Ω)
where

A : H 2 0 (Ω) → H 2 0 (Ω), B : H 2 0 (Ω) → H 2 0 (Ω), C : H 2 0 (Ω) → H 2 0 (Ω)
are defined by the Riesz representation theorem as

(5) (A u , φ) H 2 0 (Ω) := Ω 1 1 -n ∆u ∆φdx, ( 6
) (B u , φ) H 2 0 (Ω) := Ω 1 1 -n (u ∆φ + n ∆u φ)dx, and (7)
(C u , φ) H 2 0 (Ω) := Ω n 1 -n u φdx.
Note that A : H 2 0 (Ω) → H 2 0 (Ω) is a bounded, self-adjoint and invertible operator (thanks to (2)),

B : H 2 0 (Ω) → H 2 0 (Ω) is a bounded, compact and self-adjoint operator and C :

H 2 0 (Ω) → H 2 0 (Ω)
is a (non negative or non positive) bounded, compact and self-adjoint operator. Observe that since A is invertible, λ = 0. In order to avoid distinguishing the cases of 1 -n being positive or negative we shall abusively set C

1 2 ≡ -(-C 1
2 ) in the case where 1 -n non positive.

Setting U = (u , λ C 1 2 u ), the transmission eigenvalue problem (4) can be transformed into

the linear eigenvalue problem, τ ∈ R, U ∈ H 2 0 (Ω) × H 2 0 (Ω) such that (8) (T -τ I)U = 0, with τ = 1 λ ,
for the compact non-selfadjoint operator T :

H 2 0 (Ω) × H 2 0 (Ω) → H 2 0 (Ω) × H 2 0 (Ω) defined by (9) T = -A -1 B -A -1 C 1 2 C 1 2 0 .
We set (10)

T 0 = -A -1 0 B 0 -A -1 0 C 1 2 0 C 1 2 0 0
where A 0 , B 0 and C 0 are defined by ( 5), ( 6) and ( 7) respectively for n = n 0 in Ω. We state here the main result of this paper which will be proven in Section 4. In the following a transmission eigenvalue λ 0 is called simple if the corresponding τ 0 = 1/λ 0 has an algebraic multiplicity equal to 1. We refer to Theorem 4.11 for the case where λ 0 has an associated eigenspace formed only by eigenvectors (and therefore an algebraic multiplicity that coincides with the geometrical multiplicity).

Theorem 2.1. Assume that n 0 , n 1 ∈ C 4 (Ω). Let λ 0 ∈ R be a simple transmission eigenvalue of (3) with n = n 0 in Ω and let u 0 ∈ H 2 0 (Ω) be an associated eigenfunction. This implies in particular that

β 0 := Ω 1 1 -n 0 λ 2 0 n 0 |u 0 | 2 -|∆u 0 | 2 dx = 0.
If we suppose in addition that u 0 and A -1 0 u 0 are in C 6 (Ω), then, for sufficiently small > 0, there exists a transmission eigenvalue λ of (3) such that

λ = λ 0 + λ 1 + O( 3 2
)

where λ 1 is given by the following expression

λ 1 := λ 0 β 0 Γ n 0 -n 1 (1 -n 0 ) 2 |∆u 0 | 2 ds(x).
This theorem is an immediate consequence of Theorem 4.8 that is stated and proven in the last section of this paper.

The formal calculations in Section 3 show that the formula for λ 1 is generically valid whenever β 0 = 0. However, we remark that in the case of transmission eigenvalues with multiplicity greater than 1, this is not automatically ensured (See Theorem 4.11 for a rigorous expression of λ 1 that involves all eigenvectors associated with λ 0 ).

From the practical point of view, this theorem implies in particular that λ 1 gives a measure for the contrast n 0 -n 1 . For instance, if n 1 is constant and n 0 is constant on Γ, one can approximate the value of n 1 using the identity ( 11)

n 1 | Γ = n 0 | Γ - λ -λ 0 α 0 Ω 1 1 -n 0 λ 0 n 0 |u 0 | 2 - 1 λ 0 |∆u 0 | 2 dx + O( 1 2 )
with

α 0 := Γ |∆u 0 | 2 (1 -n 0 ) 2 ds(x).
For the inverse problem where one would like to determine n 1 from multistatic measurements of scattered waves, the value of λ can be approximated using sampling methods as in [START_REF] Cakoni | On the determination of dirichlet or transmission eigenvalues from farfield data[END_REF][START_REF] Cakoni | Haddar Inverse Scattering Theory and Transmission Eigenvalues SIAM publications[END_REF] (see also [START_REF] Kirsch | The inside-outside duality for scattering problems by inhomogeneous media[END_REF] for an alternative approach). The values of λ 0 and u 0 can be computed numerically if one has a priori knowledge of n 0 and Ω (see for instance [START_REF] Cakoni | Error analysis of the finite element approximation of transmission eigenvalues[END_REF][START_REF] Ji | A mixed finite element method for Helmholtz Transmission eigenvalues[END_REF][START_REF] Kleefeld | A numerical method to compute interior transmission eigenvalues Inverse Problems[END_REF] for numerical methods to approximate λ 0 and u 0 ). We finally indicate that, although not carefully checked, we conjecture that the expression for λ 1 remains true in three dimensions (corrections due to the curvature of Γ only affect higher order terms).

Formal asymptotic expansion.

In this section, we derive the formal asymptotic expansion for transmission eigenvalues and give explicit formulas for the terms up to order 2. The idea here is to provide a systematic formal way to quickly obtain the explicit expression of λ 1 in Theorem 2.1 and also higher order terms. The latter turn out to have complicated expressions that would be of marginal interest for the solution of the inverse problem mentioned above. This formal stage will also be helpful in establishing the rigorous proof based of Osborn's theorem [START_REF] Je Osborn | Spectral approximations for compacts operators[END_REF].

It allows one to have an intuition for the expression of the corrector in the asymptotic of the main operator A .

We assume the following expansions for the transmission eigenvalues :

(12) λ = ∞ j=0 j λ j ,
and then follow a classical technique for thin layers asymptotics based on rescaling and asymptotic expansion with respect to the thickness . We shall mainly follow the approach in [START_REF] Cakoni | Asymptotic analysis of transmission eigenvalue problem for a Dirichlet obstacle coated by a thin layer of a non-absorbing media[END_REF].

3.1. Scaling. We assume that the boundary Γ is C ∞ -smooth, although much less regularity is needed if we restrict ourselves to only few terms in the expansion. The issue of optimal regularity assumptions for Γ is not discussed here. However, one can check that at least a C 2 regularity is needed to get the expression of λ 1 . We parametrize Γ as

Γ = {x Γ (s), s ∈ [0, L[},
with L being the length of Γ and s is the curvilinear abscissa. At the point x Γ (s), the unit tangent vector is τ (s) := dx Γ (s) ds , the curvature κ(s) is defined by:

dτ (s) ds = -κ(s)ν(s) or dν(s) ds = κ(s)τ (s).
Within these notations, the boundary of Ω 0 is parametrized as

Γ = {x Γ (s) + ν(s), s ∈ [0, L[}.
This parametrization of the surface Γ is equivalent to the definition of Γ , for > 0 a small enough parameter.

For a function u defined in Ω , we consider ũ defined on [0, L[×]0, [ by [START_REF] Delourme | Approximate models for wave propagation across thin periodic interfaces[END_REF] ũ(s, η) = u(ϕ(s, η)) where ϕ(s, η) := x Γ (s) + ην(s).

Then, the gradient and Laplace operators are expressed in the local coordinates as:

∇u = 1 (1 + ηκ(s)) ∂ ∂s τ (s) + ∂ ∂η ν(s) ũ, (14) ∆u
= 1 (1 + ηκ) ∂ ∂s 1 (1 + ηκ) ∂ ∂s + κ (1 + ηκ) ∂ ∂η + ∂ 2 ∂η 2 ũ.
To make the formal calculations, we need to separate the thin layer and scaled it with respect to the thickness so that the equation are posed on a domain independent from . We therefore rewrite the transmission eigenvalue problem (1) in the following equivalent form (15)

                     ∆w + + k 2 n 1 w + = 0 in Ω , ∆w -+ k 2 n 0 w -= 0 in Ω 0 , ∆v + k 2 v = 0 in Ω, w + = w -, ∂w + ∂ν = ∂w - ∂ν on Γ , w + = v on Γ, ∂w + ∂ν = ∂v ∂ν on Γ.
We denote by ξ = η the stretched normal variable inside Ω and define

ϕ : G = [0, L[×]0, 1[ → Ω (s, ξ) → ϕ (s, ξ) = x Γ (s) + ξν(s).
Then the expression of the Laplace operator in the scaled layer is:

(16) ∆u = 1 (1 + ξ κ) ∂ ∂s 1 (1 + ξ κ) ∂ ∂s + κ (1 + ξ κ) ∂ ∂ξ + 1 2 ∂ 2 ∂ξ 2 û =: ∆ s,ξ û for û(s, ξ) := u(ϕ (s, ξ)).
   ŵ (s, 0) = v (x Γ (s)) s ∈ [0, L[, 1 ∂ ŵ ∂ξ (s, 0) = ∂v ∂ν (x Γ (s)) s ∈ [0, L[.
We assume that

(19) ŵ (s, ξ) = ∞ j=0 j ŵj (s, ξ), (s, ξ) ∈ G and v (x) = ∞ j=0 j v j (x), x ∈ Ω
for some functions ŵj defined on G and v j defined on Ω that are independent from . Multiplying (17) by 2 (1 + ξ κ) 3 and using [START_REF] Cakoni | Error analysis of the finite element approximation of transmission eigenvalues[END_REF], we obtain

5 k=0 k A k ŵ = 0,
where (A k ) k=0...5 are differential operators of order 2 at maximum with the following expressions for the first fourth terms:

A 0 = ∂ 2 ∂ξ 2 , A 1 =3ξκ ∂ 2 ∂ξ 2 + κ ∂ ∂ξ , A 2 = ∂ 2 ∂s 2 + 3ξ 2 κ 2 ∂ 2 ∂ξ 2 + 2ξκ 2 ∂ ∂ξ + λ 0 n 1 , A 3 =ξ 3 κ 3 ∂ 2 ∂ξ 2 + ξ 2 κ 3 ∂ ∂ξ -ξ ∂κ ∂s ∂ ∂s + ξκ ∂ 2 ∂s 2 + 3λ 0 n 1 ξκ + λ 1 n 1 .
Inserting the ansatz [START_REF] Kirsch | The inside-outside duality for scattering problems by inhomogeneous media[END_REF] in ( 17) and ( 18) we obtain after equating the terms of same order in and using the convention ŵj = v j = 0 for j < 0,

               ∂ 2 ∂ξ 2 ŵj = - 5 k=1 A k ŵj-k in G, ŵj (s, 0) = v j (x Γ (s)) s ∈ [0, L[, ∂ ŵj ∂ξ (s, 0) = ∂v j-1 ∂ν (x Γ (s)) s ∈ [0, L[. (20)
These equations can be solved inductively to get the expressions of ŵj in terms of the boundary values of v l , l ≤ j. One gets for j = 0, 1, 2 and 3

ŵ0 (s, ξ) = v 0 (x Γ (s)), ŵ1 (s, ξ) = ∂v 0 ∂ν (x Γ (s))ξ + v 1 (x Γ (s)), ŵ2 (s, ξ) = - ξ 2 2 κ ∂w - 0 ∂ν (x Γ (s)) + ∂ 2 w - 0 ∂s 2 (x Γ (s)) + λ 0 n 1 w - 0 (x Γ (s)) + ∂v 1 ∂ν (x Γ (s)ξ + v 2 (x Γ (s)), (21) and ŵ3
(s, ξ) = ξ 3 6 -2κ 2 ∂w - 0 ∂ν (x Γ (s)) -3κ ∂ 2 w - 0 ∂s 2 (x Γ (s)) -κλ 0 n 1 w - 0 (x Γ (s)) + λ 0 n 1 ∂v 0 ∂ν (x Γ (s)) + ξ 3 6 ∂ 3 v 0 ∂s 2 ∂ν (x Γ (s)) - κ ∂w - 0 ∂s (x Γ (s)) + ξ 2 2 κ ∂v 1 ∂ν (x Γ (s)) + λ 0 n 1 v 1 (x Γ (s)) + λ 1 n 1 w - 0 (x Γ (s)) + ∂v 2 ∂ν (x Γ (s))ξ + v 3 (x Γ (s)). (22)
Now, we also postulate the following expansion for w -:

(23) w -(x) = ∞ j=0 j w - j (x)
with w - j : Ω → R are functions independent of . Then (w - j , v j ) solves ( 24)

           ∆w - j + λ 0 n 0 w - j = - j l=1 λ l n 0 w - j-l in Ω, ∆v j + λ 0 v j = - j l=1 λ l v j-l in Ω.
Note that the functions w - j are defined in all Ω and not only Ω 0 and therefore [START_REF] Je Osborn | Spectral approximations for compacts operators[END_REF] gives a extension of w -to all Ω. The continuity conditions at Γ can be written as w-(s, ) = ŵ (s, 1) and

∂ w- ∂η (s, ) = 1 ∂ ŵ ∂ξ (s, 1)
where w-is defined from w -using the local change of variables ( 13) in a neighborhood of Γ.

Using Taylor's expansion (up to the second order, which is sufficient to compute the first three terms in the asymptotic expansion) we get

(25) w-(s, ) = w-(s, 0) + ∂ w- ∂η (s, 0) + 2 2 ∂ 2 w- ∂η 2 (s, 0) + o( 2 ) = ŵ (s , 1) and (26)
∂ w- ∂η (s, ) = ∂ w- ∂η (s, 0) + ∂ 2 w- ∂η 2 (s, 0) + 2 2 ∂ 3 w- ∂η 3 (s, 0) + o( 2 ) = 1 ∂ ŵ ∂ξ (s, 1).
Injecting [START_REF] Kirsch | The inside-outside duality for scattering problems by inhomogeneous media[END_REF] and ( 23) into ( 25) and ( 26), we respectively obtain the following continuity conditions on Γ, w - 0 (x Γ (s)) = ŵ0 (s, 1),

w - 1 (x Γ (s)) + ∂w - 0 ∂ν (x Γ (s)) = ŵ1 (s, 1), (27) w - 2 (x Γ (s)) + ∂w - 1 ∂ν (x Γ (s)) + 1 2 ∂ 2 w - 0 ∂ν 2 (x Γ (s)) = ŵ2 (s, 1), and 0 = ∂ ŵ0 ∂ξ (s, 1), ∂w - 0 ∂ν (x Γ (s)) = ∂ ŵ1 ∂ξ (s, 1), (28) ∂w - 1 ∂ν (x Γ (s)) + ∂ 2 w - 0 ∂ν 2 (x Γ (s)) = ∂ ŵ2 ∂ξ (s, 1).
System [START_REF] Reed | Functional analysis[END_REF] coupled with the boundary conditions (28) and ( 27) provide an inductive way to determine (w - j , v j ). We obtain the set of equations satisfied by these terms after substituting the expressions of ŵj (s, 1) given by ( 21), [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]. We hereafter summarize the set of equations obtained for (w - j , v j ) and how to use them to get the expressions of λ j , j = 0, 1, 2.

We first obtain that the couple (w - 0 , v 0 ) solves (29)

         ∆w - 0 + λ 0 n 0 w - 0 = 0 in Ω, ∆v 0 + λ 0 v 0 = 0 in Ω, w - 0 -v 0 = 0 on Γ, ∂w - 0 ∂ν - ∂v 0 ∂ν = 0 on Γ.
This means in particular that λ 0 is a transmission eigenvalue for the limiting problem where the thin layer is removed. We then obtain that the couple (w - 1 , v 1 ) satisfies (30)

         ∆w - 1 + λ 0 n 0 w - 1 = -λ 1 n 0 w - 0 in Ω, ∆v 1 + λ 0 v 1 = -λ 1 v 0 in Ω, w - 1 -v 1 = 0 on Γ, ∂w - 1 ∂ν - ∂v 1 ∂ν = λ 0 (n 0 -n 1 )w - 0 on Γ.
Since λ 0 is an eigenvalue of the associated homogeneous system, this problem is solvable only if a compatibility condition is satisfied by the right hand sides. This compatibility condition can be obtained by multiplying the first equation with w - 0 and the second equation with v 0 , taking the difference then integrating by parts and using (29). One ends up with

λ 1 = Γ λ 0 (n 0 -n 1 )|w - 0 | 2 ds(x) Ω n 0 |w - 0 | 2 -|v 0 | 2 dx
which coincides with the expression of given in Theorem 2.1 expressed in terms of u 0 = w - 0 -v 0 .

Although not covered by the analysis of convergence, we also provide the expression of the third term in the asymptotic expression. One get that the couple (w - 2 , v 2 ) solves (31)

         ∆w - 2 + λ 0 n 0 w - 2 = -λ 1 n 0 w - 1 -λ 2 n 0 w - 0 in Ω, ∆v 2 + λ 0 v 2 = -λ 1 v 1 -λ 2 v 0 in Ω, w - 2 -v 2 = h 1 on Γ, ∂w - 2 ∂ν - ∂v 2 ∂ν = h 2 on Γ,
where

h 1 = - 1 2 ∂ 2 w - 0 ∂ν 2 - 1 2 λ 0 (n 0 -n 1 )w - 0 and h 2 = κ ∂ 2 w - 0 ∂ν 2 - 7κ 2
∂ 2 w - 0 ∂s 2 + 2κ 2 + λ 0 (n 0 + n 1 2 ) ∂w - 0 ∂ν - 3κ 2
∂w - 0 ∂s + 3 2 ∂ 3 w - 0 ∂ν∂s 2 + λ 1 (2n 1 -n 0 ) + λ 0 (κ( n 1 2 -n 0 )) w - 0 - ∂ 2 w - 1 ∂ν 2 + κ ∂w - 1 ∂ν .
Writing the compatibility condition for (31), we obtain the following formula for λ 2

λ 2 Ω 1 1 -n 0 1 λ 0 |∆u 0 | 2 -λ 0 n 0 |u 0 | 2 dx = -λ 2 1 Ω 1 λ 0 ∆u 0 ū0 + 1 1 -n 0 |u 0 | 2 dx -λ 1 Ω 1 1 -n 0 u 1 ∆ū 0 + n 0 ∆u 1 ū0 + 2n 0 λ 0 u 1 ū0 dx + Γ h 1 ∂ ∂ν 1 1 -n 0 (∆ + λ 0 )ū 0 ds(x) - Γ h 2 1 1 -n 0 (∆ + λ 0 )ū 0 ds(x). ( 32
)
This complicated expression shows in particular a nonlinear dependence of λ 2 in terms of n 1 . It suggests that the use of λ 2 for solutions to the inverse problem of determining n 1 may not be appropriate.

4. Convergence analysis. The main goal of this section is to prove Theorem 2.1 that provides a rigorous mathematical justification of the formal asymptotic expansion for simple real transmission eigenvalues up to the first order. The proof is split into several steps. The first one is to establish the convergence in norm of the operator T to T 0 . This ensures the convergence of λ to λ 0 . In order to get to the term of order 1 in , we shall apply the Osborn theorem which requires for instance a characterization of the pointwise asymptotic expansion of T (U ) up to order 1 in (for some given function U ∈ H 2 0 (Ω) × H 2 0 (Ω)). The latter can be obtained from the asymptotic expansions of A -1 u, B u and C u for some u ∈ H 2 0 (Ω). The difficult part to get the expansion of A -1 u since for the two others, the first order terms are vanishing. This critical result is provided by Lemma 4.5.

In all the following we use the notation

(f, g) := (f, g) H 2 0 (Ω) = Ω
∆f ∆gdx and g := (g, g)

1 2
H 2 0 (Ω) .

For an operator A : V → V , A denotes the operator norm. To simplify the writing, C will denote a generic constant whose value may change but remains independent from as → 0.

4.1. Pointwise convergence of the spectrum of T . In this first step, we prove pointwise convergence of the spectrum of the operator T to the spectrum of T 0 . This is a direct consequence of the following convergence in norm [START_REF] Je Osborn | Spectral approximations for compacts operators[END_REF][START_REF] Cakoni | Asymptotic expansions for transmission eigenvalues for media with small inhomogeneities[END_REF].

Theorem 4.1. Assume that n 0 ∈ C 2 (Ω). Let T and T 0 be defined by ( 9) and (10) respectively.

Then T converges to T 0 in the operator norm.

Proof. The proof follows from Lemma 4.2 and Lemma 4.4 below, using the definition of T and T 0 .

In the first lemma we prove norm convergence for B and C .

Lemma 4.2. Let B , C , B 0 and C 0 be the operators defined by ( 6) and [START_REF] Cakoni | The existence of an infinite discrete set of transmission eigenvalues[END_REF]. Then, for sufficiently small ,

(33) B -B 0 ≤ C 1 2
and C

1 2 -C 1 2 0 ≤ C .
Proof. From the definitions of B and B 0 , we have that for u, φ ∈ H 2 0 (Ω)

((B -B 0 )u, φ) = Ω 1 1 -n u∆φ + n ∆uφ dx - Ω 1 1 -n 0 u∆φ + n 0 ∆uφ dx = Ω 1 1 -n 1 - 1 1 -n 0 u∆φ + ∆uφ dx.
Therefore,

|((B -B 0 )u, φ)| ≤C u L ∞ (Ω) ∆φ L 1 (Ω ) + φ L ∞ (Ω ∆u L 1 (Ω ) .
Using the Sobolev embedding theorem and the Cauchy Schwartz inequality, we get

|((B -B 0 )u, φ)| ≤ C 1 2 ( u H 2 0 (Ω) φ H 2 0 (Ω) ).
By choosing φ = (B -B 0 )u, we get

(B -B 0 )u H 2 0 (Ω) ≤ C 1 2 u H 2 0 (Ω) .
The proof is similar for the second inequality. For u, φ ∈ H 2 0 (Ω), we have

((C -C 0 )u, φ) = Ω 1 1 -n 1 - 1 1 -n 0 uφdx ≤ C |Ω | u L ∞ (Ω) φ L ∞ (Ω)
From the Sobolev embedding theorem, we obtain

((C -C 0 )u, φ) ≤ C u H 2 0 (Ω) φ H 2 0 (Ω) .
By choosing φ = (C -C 0 )u, we have

(34) (C -C 0 )u H 2 0 (Ω) ≤ C u H 2 0 (Ω) .
Using the square root Lemma in [START_REF] Reed | Functional analysis[END_REF] and the fact that C n converges to C n 0 at the same order O( ), we conclude that C

1 2 converges to C 1 2
0 at the same order O( ). Hence we have

(35) (C 1 2 -C 1 2 0 )u H 2 0 (Ω) ≤ C u H 2 0 (Ω) .
Now we show convergence in the H 2 0 (Ω) norm for A -1 f assuming smoothness of f . This will be useful in the proof of Lemma 4.4 since the operators B and C are regularizing.

Lemma 4.3. Let A and A 0 be defined by (5) for > 0 and = 0, respectively and f ∈ H 2 0 (Ω).

If A -1 0 f ∈ C 2 (Ω), then for sufficiently small , (36)

A -1 f -A -1 0 f ≤ C 1 2 .
Proof. For a fixed f ∈ H 2 0 (Ω), define z and z 0 in H 2 0 (Ω) as z = A -1 f and z 0 = A -1 0 f . Since

A z = A 0 z 0 = f , we have that for φ ∈ H 2 0 (Ω) (37) (A (z -z 0 ), φ) = (A 0 z 0 -A z 0 , φ) = Ω 1 1 -n 0 - 1 1 -n 1 ∆z 0 ∆φdx. If z 0 ∈ C 2 (Ω), we get Ω 1 1 -n 0 - 1 1 -n ∆z 0 ∆φdx ≤ C ∆z 0 ∞ Ω ∆φdx ≤C 1 2 φ H 2 0 (Ω) .
Thus, we have shown that

(A (z -z 0 ), φ) ≤ C 1 2 φ H 2 0 (Ω) .
By plugging in φ = z -z 0 , we obtain the desired convergence using the coercivity of A .

Lemma 4.4. Assume that n 0 ∈ C 2 (Ω). Let A , B , C , A 0 , B 0 and C 0 be defined by ( 5), ( 6) and (7) for > 0 and = 0, respectively. Then for sufficiently small ,

A -1 B -A -1 0 B 0 -→ →0 0 and A -1 C 1 2 -A -1 0 C 1 2 0 -→ →0 0.
Proof. From (37), we have that for f, φ ∈ H 2 0 (Ω) and with z = A -1 f and

z 0 = A -1 0 f (A (z -z 0 ), φ) ≤ C ∆A -1 0 f L 2 (Ω ) φ H 2 0 (Ω) .
Furthermore,

A -1 B f -A -1 0 B 0 f H 2 0 (Ω) ≤ (A -1 -A -1 0 )B 0 f H 2 0 (Ω) + A -1 (B -B 0 )f H 2 0 (Ω) ≤C ∆A -1 0 B 0 f L 2 (Ω ) + A -1 (B -B 0 ) f H 2 0 (Ω) . (38) For estimating ∆A -1 0 B 0 f L 2 (Ω ) , observe that B 0 u ∈ H 2 0 (Ω) is the weak solution ∆∆B 0 u = ∆ n 0 1 -n 0 u + 1 1 -n 0 ∆u in Ω.
Classical regularity results [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF][START_REF] Wloka | Partial Differential Equations[END_REF] and the fact that

n 0 ∈ C 2 (Ω) imply that B 0 u ∈ H 4 (Ω) ∩ H 2 0 (Ω)
and therefore

∆A -1 0 B 0 f H 1 (Ω) ≤ C f H 2 (Ω) .
By the Sobolev embedding theorem, this implies that

∆A -1 0 B 0 f L p (Ω) ≤ C f H 2 (Ω) ,
for p > 2 . Let p = p 2 > 1 and q such that

1 p + 1 q = 1. ∆A -1 0 B 0 f 2 L 2 (Ω ) ≤ ∆A -1 0 B 0 f 2 L p (Ω) |Ω | 1 q ≤ C 1 2q f H 2 (Ω) . (39)
From (33) we obtain (40)

A -1 (B -B 0 ) f H 2 0 (Ω) ≤ C 1 2 f H 2 0 (Ω) .
Using (38), ( 39) and (40) we have that

A -1 B -A -1 0 B 0 -→ →0 0.
The second convergence result follows from similar arguments. Now we would like to obtain explicit formula for the correction term in the asymptotic expansion for the operator T . More precisely, we define explicit formula for the corrector term associated with A -1 -A -1 0 .

Corrector term for

A -1 -A -1 0 .
In this subsection, we construct a corrector function and use it to estimate the convergence rate of z

= A -1 u for u ∈ H 2 0 (Ω). Let z 0 = A -1 0 u ∈ H 2 0 (Ω), i.e z 0 ∈ H 2 0 (Ω) solution of (41) ∆ 1 1 -n 0 ∆z 0 = ∆∆u in Ω.
Inspired by the formal calculations on the previous section, namely problem (30), we define z 1 solution of (42)

         ∆ 1 1 -n 0 ∆z 1 = 0 in Ω, z 1 = 0 on Γ, ∂z 1 ∂ν = 1 -n 1 1 -n 0 -1 ∆z 0 on Γ.
We expect that z = z 0 + z 1 + O( 2) in Ω 0 . We extend z 1 in Ω as z 1 defined by ( 43)

z 1 = z 1 in Ω 0 , z 1 -ψ in Ω
where ψ is a polynomial of order ≤ 3 and satisfying the boundary conditions:

(44)      ψ = 0 , ∂ψ ∂ν = 1 -n 1 1 -n 0 -1 ∆z 0 on Γ, ψ = ∂ψ ∂ν = 0 on Γ .
This gives the following expression of ψ (that plays the role ŵ2 in the formal calculations)

ψ(x) = ψ(ϕ(s, ξ)) = ψ(s, ξ) = 1 -n 1 1 -n 0 -1 ∆z 0 (ϕ(s, 0))ξ(1 -ξ) 2 .
The choice of ψ ensures in particular that z 1 ∈ H 2 0 (Ω). To simplify the notation we set

m := 1 1 -n 0 - 1 1 -n 1 .
Now we have the following Lemma.

Lemma 4.5. Assume that n 0 and n 1 are in C 4 (Ω). Let u ∈ H 2 0 (Ω) then set z = A -1 u and z 0 = A -1 0 u. We define z 1 as in (43) and assume that z 0 ∈ C 6 (Ω). Then we have, for sufficiently small ,

z -z 0 -z 1 H 2 0 (Ω) ≤ C 3 2 .
Proof. For any φ ∈ H 2 0 (Ω) we have that

(45) (A (z -z 0 -z 1 ), φ) = (A (z -z 0 ), φ) -(A z 1 , φ).
We recall that

(A (z -z 0 ), φ) = Ω 1 1 -n 0 - 1 1 -n 1 ∆z 0 ∆φdx.
Furthermore, we have that

(46) (A z 1 , φ) = Ω 0 1 1 -n 0 ∆z 1 ∆φdx + Ω 1 1 -n 1 ∆(z 1 -ψ)∆φdx.
Using the fact that ∆ 1 1 -n 0 ∆z 1 = 0 and the Green formula yield,

(A z 1 , φ) = Γ m∆z 1 ∂φ ∂ν ds(x) + Γ ∂ ∂ν 1 1 -n 1 ∆z 1 - ∂ ∂ν 1 1 -n 0 ∆z 1 φds(x) + Γ 1 1 -n 1 ∆ψ ∂φ ∂ν ds(x) - Γ ∂ ∂ν 1 1 -n 1 ∆ψ φds(x) - Ω ∆ 1 1 -n 1 ∆ψφdx.
Using the expression of ψ we have

1 1 -n 1 ∆ψ| Γ = 1 1 -n 1 ∆ s,ξ ψ(s, 1) = 2 m∆z 0 (ϕ(s, 0)), ∂ ∂ν ( 1 1 -n 1 ∆ψ)| Γ = ∂ ∂η ( 1 1 -n 1 ∆ψ)| Γ = 1 1 -n 1 1 ∂ ∂ξ (∆ s,ξ ψ)(s, 1) = 6 2 m∆z 0 (ϕ(s, 0)).
We then get after substitution of these expressions

(A z 1 , φ) = Γ m ∆z 1 (ϕ(s, )) + 2 ∆z 0 (ϕ(s, 0)) ∂φ ∂ν ds(x) - Γ m ∂ ∂ν (∆z 1 (ϕ(s, ) + 6 2 ∆z 0 (ϕ(s, 0)))φds(x) - Ω ∆ 1 1 -n 1 ∆ψφdx = Γ φ 1 ∂φ ∂ν ds(x) - Γ φ 2 φds(x) - Ω ∆ 1 1 -n 1 ∆ψφdx
where we have set

φ 1 (s) := m∆z 1 (ϕ(s, )) + 2 m∆z 0 (ϕ(s, 0)), φ 2 (s) := m ∂ ∂ν (∆z 1 (ϕ(s, ))) + 6 2 m∆z 0 (ϕ(s, 0))
using the parametrization of the curve Γ , s -→ ϕ(s, ) with ϕ defined by [START_REF] Delourme | Approximate models for wave propagation across thin periodic interfaces[END_REF]. Using this parametrization and setting φ(s, η) := φ(ϕ(s, η)) in Ω we have

Γ φ 1 ∂φ ∂ν ds(x) = L 0 φ 1 ∂ φ ∂η (s, )(1 + κ)ds = L 0 0 φ 1 ∂ 2 φ ∂η 2 (s, η)(1 + κ)dsdη.
From the definition of φ 1 we then get for φ ∈ H 2 0 (Ω),

Γ φ 1 ∂φ ∂ν ds(x) = 2 L 0 0 m∆z 0 (ϕ(s, 0)) ∂ 2 φ ∂η 2 (s, η)dsdη + O( 1
2 ) φ H 2 (Ω) .
Here and in all the following O( r ) denotes a function such that O( r ) ≤ C r for a constant C independent from the test function φ but that may depend on z 0 C 6 (Ω) . Using Taylor's expansion we also get for φ ∈ H 2 0 (Ω),

Γ φ 2 φds(x) = 2 L 0 0 φ 2 ∂ 2 φ ∂η 2 (s, η)(1 + κ)dsdη + O( 3
2 ) φ H 2 (Ω) = 3 L 0 0 m∆z 0 (ϕ(s, 0)) ∂ 2 φ ∂η 2 (s, η)dsdη + O( 1 2 ) φ H 2 (Ω)
where the last equality is obtained after substituting the expression of φ 2 . One ends up with

Γ (φ 1 ∂φ ∂ν -φ 2 φ)ds(x) = - L 0 0 m∆z 0 (ϕ(s, 0)) ∂ 2 φ ∂η 2 (s, η)dsdη + O( 3
2 ) φ H 2 (Ω) .
Equation (45) then gives

(A (z -z 0 -z 1 ), φ) = Ω m∆z 0 ∆φdx -(A z 1 , φ) = L 0 0 m∆z 0 (ϕ(s, η))∆φ(ϕ(s, η))(1 + ηκ)dsdη - L 0 0 m∆z 0 (ϕ(s, 0) ∂ 2 φ ∂η 2 (ϕ(s, η))dsdη - Ω ∆ 1 1 -n 1 ∆ψφdx + O( 3
2 ) φ H 2 (Ω) .
We use the expression of the Laplacien in local coordinates

(1 + ηκ)∆φ(ϕ(s, η)) = ∂ ∂s 1 1 + ηκ ∂ φ ∂s (s, η) + κ ∂ φ ∂η (s, η) + (1 + ηκ) ∂ 2 φ ∂η 2 (s, η)
to make the decomposition

L 0 0 m∆z 0 (ϕ(s, η))∆φ(ϕ(s, η))(1 + ηκ)dsdη = L 0 0 m∆z 0 (ϕ(s, η)) ∂ ∂s 1 1 + ηκ ∂ φ ∂s (s, η) + L 0 0 m∆z 0 (ϕ(s, η)) κ ∂ φ ∂η (s, η) + ηκ ∂ 2 φ ∂η 2 (s, η) + L 0 0 m∆z 0 (ϕ(s, η)) ∂ 2 φ ∂η 2 (s, η).
To estimate the first term, we integrate by parts on [0, L[, we obtain

L 0 0 m∆z 0 (ϕ(s, η)) ∂ ∂s 1 1 + ηκ ∂ φ ∂s (s, η) dηds = - L 0 0 1 1 + ηκ ∂ ∂s (m∆z 0 (ϕ(s, η))) ∂ φ ∂s (s, η)dsdη = - L 0 1 0 m 1 1 + ξκ ∂ ∂s ∆z 0 (ϕ(s, ξ)) ∂ φ ∂s (s, ξ)dsdξ = - L 0 1 0 m ∂ ∂s ∆z 0 (ϕ(s, 0)) ξ 0 ∂ 2 φ ∂η∂s (s, η)dη dsdξ + O( 2 ) φ H 2 (Ω) = O( 3 2
) φ H 2 (Ω) .

For the last term we proceed similarly to obtain

L 0 0 m∆z 0 (ϕ(s, η)) ∂ φ ∂η (s, η)dsdη = L 0 1 0 m∆z 0 (ϕ(s, ξ)) ∂ φ ∂η (s, ξ) dsdξ = L 0 1 0 m∆z 0 (ϕ(s, 0)) ξ 0 ∂ 2 φ ∂η 2 (s, η)dη dsdξ + O( 2 ) φ H 2 (Ω) = O( 3
2 ) φ H 2 (Ω)
Observing in addition that

L 0 0 ηκm(∆z 0 (ϕ(s, η) ∂ 2 φ ∂η 2 (s, η))dsdη = O( 3 2
) φ H 2 (Ω) , one ends up with

(A (z -z 0 -z 1 ), φ) = L 0 0 m ∆z 0 (ϕ(s, η)) -∆z 0 (ϕ(s, 0)) ∂ 2 φ ∂η 2 (s, η)dsdη - Ω ∆ 1 1 -n 1 ∆f ψ φdx + O( 3 2
) φ H 2 (Ω) .

To conclude we just observe that the two remaining terms are also of the form O(

) φ H 2 (Ω) . For the first term, we simply use a Taylor expansion for ∆z 0 while for the second one we just use that, due to the regularity of n 0 and n 1 ,

∆ 1 1 -n 1 ∆ψ ∈ L ∞ (Ω).
In conclusion,

(A (z -z 0 -z 1 ), φ) ≤ C 3 2 φ H 2 (Ω) .
Choosing φ = z -z 0 -z 1 , since the coercivity constant associated with A is independent from , we get

z -z 0 -z 1 H 2 0 (Ω) ≤ C 3 2
which ends the proof.

B 0 (A -1 -A -1 0 )u H 2 0 (Ω) ≤ C and C 1 2 0 (A -1 -A -1 0 )u H 2 0 (Ω) ≤ C (47)
where C independent of .

Using the local coordinates in Ω , we obtain

Ω 1 1 -n 1 - 1 1 -n 0 ∆uφdx = L 0 1 0 m∆u(ϕ(s, ξ)) φ(s, ξ) (1 + ξ κ)dsdξ = L 0 1 0 m∆u(ϕ(s, 0)) ξ 0 ∂ 2 φ ∂η 2 (s, η)dη (1 + ξ κ)dsdξ ≤ C 3 2 φ H 2 (Ω) Hence ((B -B 0 )u, φ) ≤ C 3 2 φ H 2 (Ω)
. Similarly, we compute the asymptotic formula of

((C -C 0 )u, φ) = L 0 1 0 n 0 1 -n 0 u(ϕ(s, ξ))φ(ϕ(s, ξ)) (1 + ξ κ)dsdξ = L 0 1 0 n 0 1 -n 0 ξ 0 ∂ ũ ∂η (s, η)dη ξ 0 ∂ φ ∂η (s, η)dη (1 + ξ κ)dsdξ ≤ C 2 φ H 2 (Ω)
Using the square root Lemma in [START_REF] Reed | Functional analysis[END_REF] and the fact that C n converges to C n 0 at the same order

O( 2 ), we conclude that C 1 2 converges to C 1 2
0 at the same order O( 2). Thus we have

((C 1 2 -C 1 2 0 )u, φ) ≤ C 2 φ H 2 (Ω)
Consequently

and concludes the proof.

We now indicate a possible extension to the case where the eigenvalue τ 0 is not simple. We need in that case to assume that the geometric multiplicity m coincides with the algebraic multiplicity so that a basis of R(E) is formed by eigenvectors of T 0 that we denote by U j 0 = u j 0 λ 0 C . If we assume that (60) -β j 0 := U j 0 , Ũ j * 0 = (A 0 u j 0 , u j 0 ) -λ 2 0 (C 0 u j 0 , u j 0 ) = 0, then we can define the dual basis as U j * 0 = -1

β j 0 Ũ j * 0 .
Notice that the assumption on β j 0 is not guaranteed in general. Making this assumption makes the expression of the dual basis easier to express and allows us to follow the same calculations as above to express the leading term in (T 0 -T )U j 0 , U j * 0 .

We then obtain the following result as a consequence of the application of Theorem 4.7.

Theorem 4.11. Assume that n 0 and n 1 are in C 4 (Ω). Let λ 0 be a real transmission eigenvalue corresponding to n 0 such that the associated eigenspace is formed only with eigenvectors u j 0 ∈ H 2 0 (Ω), j = 1, . . . m. Assume in addition that β j 0 defined by (60) does not vanish and that u j 0 and A -1 0 u j 0 are in C 6 (Ω). Then, for > 0 small enough, there exists m transmission eigenvalues λ j corresponding to n such that

(61) 1 m m i=1 1 λ j - 1 λ 0 = - λ 0 1 m m i=1 1 β j 0 Γ n 0 -n 1 (1 -n 0 ) 2 |∆u j 0 | 2 ds(x) + O( 3 2
).

Fig. 1 .

 1 Fig. 1. Stretch of the geometry

  The next step is to write the equation for w + in the scaled domain and solve for the asymptotic expansion of w + in terms of the boundary values on Γ. These boundary values are given by the asymptotic expansion of v . More specifically, setting ŵ (s, ξ) := w + (ϕ (s, ξ)), we have that[START_REF] Henrot | Extremum Problems for Eigenvalues of Elliptic Operators[END_REF] ∆ s,ξ ŵ + λ n 1 ŵ = 0 in G together with the boundary conditions[START_REF] Ji | A mixed finite element method for Helmholtz Transmission eigenvalues[END_REF]

Lemma 4 . 6 .

 46 Assume that n 0 and n 1 are in C 4 (Ω). If u ∈ C 6 (Ω) ∩ H 2 0 (Ω), then for sufficiently small ,

Lemma 4 . 10 . 0 n 0 -n 1 ( 1 - 2 ) 1 1 0 m 1 -n 1 1 0 m 1 -n 1 1

 41001121011011 Under the assumptions of Theorem 4.8 one has for any φ ∈ H2 0 (Ω) ∩ C 4 (Ω) (A 0 (A -1 -A -1 0 )u, φ) = L n 0 ) 2 ∆z 0 (ϕ(s, 0))∆φ(ϕ(s, 0))ds + O( 3 where z 0 := A -1 0 u 0 . Proof. With z := A -1 u 0 , -z 0 -z 1 )∆φdx + Ω m∆(z 0 + z 1 )∆φdxApplying Lemma 4.5 we obtainΩ -n 0 ∆(z -z 0 )∆φdx = Ω m∆(z 0 + z 1 )∆φdx + O( ϕ(s, ξ))∆φ(ϕ(s, ξ)) (1 + ξκ)dsdξ = L 0 m∆z 0 (ϕ(s, 0))∆φ(ϕ(s, 0))ds + O( ξ)∆m∆φ(ϕ(s, ξ)) (1 + ξκ)dsdξ -L -n 0 -1∆z 0 (ϕ(s, 0))∆φ(ϕ(s, ξ))ds = L -n 0 -1 ∆z 0 (ϕ(s, 0))∆φ(ϕ(s, 0))ds + O(

Ω 1 1 n 1 ( 1 -

 111 -n 0 ∆(z -z 0 )∆φdx = L 0 n 0 -n 1 (1 -n 0 ) 2 ∆z 0 (ϕ(s, 0))∆φ(ϕ(s, 0))ds + O( n 0 ) 2 ∆z 0 (ϕ(s, 0))∆φ(ϕ(s, 0))ds + O(

).

Proof. From the estimate of Lemma 4.5 we have that

Since B 0 is two orders smoothing, we have that

The same proof holds for C 0 .

Now to derive the eigenvalue expansion, we will apply the Theorem of Osborn [START_REF] Je Osborn | Spectral approximations for compacts operators[END_REF], which we state here for reader's convenience. Suppose X is a Hilbert space with inner product denoted by , and K n : X → X is a sequence of compact linear operators such that K n converge in the operator norm to K. It then follows that the adjoint operators also converges in norm. Let µ be a nonzero eigenvalue of K of algebraic multiplicity m . It is well known that for n large enough, there exist m

Let E be the spectral projection onto the generalized eigenspace of K corresponding to eigenvalue µ. The space X can be decomposed in terms of the range and null space of E as X = R(E) ⊕ N (E). Then form the proof of Theorem 3 in [START_REF] Je Osborn | Spectral approximations for compacts operators[END_REF], one can state the following theorem.

Theorem 4.7. Let φ 1 , φ 2 , ..., φ m be a normalized basis for R(E), and let φ * 1 , φ * 2 , ..., φ * m be the dual basis of R(E) such that v, φ * j = 0 for all v ∈ N (E). Then there exists a constant C such that :

In order to apply Theorem 4.7 and obtain explicit expression for the first order asymptotic m j=1 (K -K n )φ j , φ * j , one has to construct the basis φ * j . Remark that in the case of selfadjoint operators, φ * j = φ j , but this does not apply to our problem. One easily check that φ * j are necessarily a basis of the generalized eigenspace of K * associated with µ.

We now turn our attention to application of this theorem with K n ≡ T and K ≡ T 0 and X ≡ H 2 0 (Ω) × H 2 0 (Ω). We already showed that T converges to T 0 in the operator norm in Theorem 4.1. In order to simplify the calculations we define the inner product on H 2 0 (Ω) × H 2 0 (Ω) by:

Let τ 0 be a simple real eigenvalue of T 0 , then for small enough, some eigenvalue

0 u 0 be an eigenvector of T 0 associated with τ 0 . Using the expression of T 0 one

0 u 0 is an eigenvector of T * 0 associated with τ 0 . Then this eigenvector is proportional to the dual basis of U 0 if and only if

We remark that since τ 0 is assumed to be a simple eigenvalue (i.e. also with geometrical multiplicity equals 1), then (49) holds. We then can define the dual vector as

and apply Theorem 4.7 to get that (50)

We are now in position to prove the main result of Theorem of this paper. We refer to Theorem 4.11 for an extension to the case of transmission eigenvalues with higher multiplicities.

Theorem 4.8. Assume that n 0 and n 1 are in C 4 (Ω). Let λ 0 be a simple real transmission eigenvalue corresponding to n 0 and let u 0 ∈ H 2 0 (Ω) be the corresponding eigenvector. This implies in particular that (49) holds. Further assume that u 0 and A -1 0 u 0 are in C 6 (Ω). Then, for > 0 small enough, there exists a transmission eigenvalue λ corresponding to n such that

).

Proof. Using estimate (50) with λ 0 = 1 τ 0 , we have

From the definition of ( 9) of T , we have

Using the definition (10) of T 0 , we obtain

0 u 0 0 From (33) and (36), we have 2 ).

On the other hand, we have

From estimates (33) and (47), we obtain

Next, (52) implies

Using the expression of U * 0 we see that

Recall that, by definition of u 0 (55)

0 is self-adjoint, we have

).

We then deduce

).

In order to conclude, we use the results of the two lemmas below that treat the asymptotic for each term in (56).

Applying Lemma 4.9 with u = u 0 and φ = u 0 we have

and with u = u 0 and φ = C

Applying Lemma 4.10 with u = A 0 u 0 (therefore z 0 = u 0 ) and φ = u 0 , we get, using the fact that A 0 is selfadjoint,

(1 -n 0 ) 2 ∆u 0 (ϕ(s, 0))∆u 0 (ϕ(s, 0))ds + O(

We finally obtain for some C independent of and φ.

Proof. Since