Yasmina Fedala

Sorin Munteanu

Frédéric Kanoufi

Gilles Tessier

Jean Paul Roger

Chang Wu

Fabien Amiot
email: fabien.amiot@femto-st.fr

Calibration procedures for quantitative multiple wavelengths reflectance microscopy

In order to characterize surface chemo-mechanical phenomena driving micro-electromechanical systems (MEMS) behavior, it has been previously proposed to use reflected intensity fields obtained from a standard microscope for different illumination wavelengths. Wavelength-dependent and -independent reflectivity fields are obtained from these images provided the relative reflectance sensitivities ratio can be identified. This contribution focuses on the necessary calibration procedures and mathematical methods allowing for a quantitative conversion from a mechanically-induced reflectivity field to a surface rotation field, therefore paving the way for a quantitative mechanical analysis of MEMS under chemical loading.

I. INTRODUCTION

Because of their high surface over volume ratio, the surface of micrometer sized structures plays a key role in their mechanical behavior, and this property has been proposed to devise micromechanical sensors of environmental changes 1 . A significant effort has been put on the development of biological sensors 2 , but the basic understanding of coupled surface phenomena [START_REF] Hagan | Nanomechanical Forces Generated by Surface Grafted DNA[END_REF] has somehow been left behind.

Imaging techniques usually reveal the heterogeneity of surface affinities at the micrometer scale [START_REF] Munteanu | Mapping fluxes of radicals from the combination of electrochemical activation and optical microscopy[END_REF] , so that the assessment of the chemical homogeneity of the considered surface is required when focusing on chemo-mechanical coupling phenomena. As surface chemical composition modifications induce (complex) reflection coefficient changes, Jin et al. [START_REF] Jin | A biosensor concept based on imaging ellipsometry for visualization of biomolecular interactions[END_REF] proposed an ellipsometric imaging set-up to measure the optical thickness of thin adsorbed films, when Li et al. [START_REF] Li | Interferometric Measurement of Depletion Layer Structure and Voltammetric Data in Concentrated Organic Redox Solutions[END_REF] use interferometry to measure locally the concentration profiles of reactants near an electrode. Differential reflectance changes detections have been achieved to perform sensitive electroreflectance [START_REF] Feinleib | Electroreflectance in metals[END_REF] and thermoreflectance 8 measurements, providing access to electrochemically induced effects at substrate-electrolyte interfaces and surface temperature fields.

Focusing on chemo-mechanical coupled phenomena also requires to measure kinematic fields of deformable surfaces. In the field of cantilever sensors, Mertens et al. [START_REF] Mertens | Real-time profile of microcantilevers for sensing applications[END_REF] coupled a scanning stage to the standard optical lever technique to measure cantilever sensors profiles.

Several interferometric techniques have also been proposed relying on phase changes related to the out-of-plane displacement field [START_REF] Amiot | Nomarski imaging interferometry to measure the displacement field of MEMS[END_REF] . This optical lever technique could also be combined to ellipsometric measurements in order to simultaneously monitor the mean surface curvature and average molecules surface densities on the surface of a microcantilever [START_REF] Godin | Combined in situ Micromechanical Cantilever-Based Sensing and Ellipsometry[END_REF] .

In order to retrieve spatially resolved data, we earlier proposed to use multiple wavelengths imaging reflection microscopy [START_REF] Garraud | Multiple wavelengths reflectance microscopy to study the multi-physical behavior of MEMS[END_REF] . This technique only requires minor modifications of commercially available reflection microscopes and can provide both local surface modification and kinematic field measurements. A decoupling method has been presented to distinguish wavelength-dependent (e.g., surface changes) and -independent (e.g., kinematic) contributions to the collected intensity. This decoupling method requires the identification of the ratio of the relative reflectance sensitivities for the considered wavelengths. In addition, converting the wavelength-independent contributions into a quantitative surface deformation requires a dedicated calibration procedure. This paper therefore aims at providing the details of robust but simple and straightforward procedures allowing one to retrieve the calibration parameters necessary to decouple wavelength-dependent and wavelength-independent contributions and to provide quantitative kinematic data. The set-up as well as the decoupling procedure are first detailed. The procedure for the identification of the relative reflectance sensitivities ratio is described and simple criteria to assess the identification procedure quality are proposed. The calibration method used to convert the wavelength-independent reflectivity change field into mechanical information is then detailed, and the conversion method is also presented, possibly taking chromatic aberrations into account. Results obtained when investigating the electro-elastic coupling on a microcantilever beam immersed in an electrolyte are used to illustrate the different steps throughout the paper.

II. MULTIPLE WAVELENGTHS REFLECTANCE MICROSCOPY

A. Samples

The deformable mechanical microcantilevers used throughout the paper to exemplify the different calibration and measurement procedures are made of LPCVD silicon nitride whose stoichiometry is adjusted to minimize residual stresses. The considered microcantilever is 431 nm thick, 80 micrometers long and 15 micrometers wide. This microcantilever (see Fig. 1) is covered with a 10 nm thick chromium adhesion layer and a gold electrode (50 nm) [START_REF] Wu | Design and fabrication of a multiple-thickness electrochemical cantilever sensor[END_REF] . The electrodes are such that each cantilever is electrically addressed independently, and a SU8 insulation layer is used to ensure that only the cantilever surface in addition to a 15 × 12µm 2 electrode area on the substrate are in constant contact with the surrounding solution. The temperature of the cantilever is controlled within 0.1 • C (closed-loop temperature regulation of a Peltier plate) for all the experimental results reported herein.

The described procedures apply to any phenomena modifying the surface reflectivity and this paper presents results obtained when involving the electro-elastic coupling. These are obtained by placing a degassed KCl 1×10 -2 M solution around the cantilever and by controlling the electrode potential (working electrode) with a home-made electrochemical workstation [START_REF] Gopinath | An inexpensive field-portable programmable potentiostat[END_REF] . A platinum wire is used as a counter-electrode, and another one is used as a pseudo-reference electrode.

B. Set-up

The objects are observed with an objective lens (Olympus UMPLFLN ×20, NA = 0.5) and imaged on a CCD array (Dalsa 1M30, 12 bits, 1024 × 1024 pixels) using focusing optics (focal length 200 mm) (Fig. 2). In order to distinguish the wavelength-dependent andindependent contributions, the sample is illuminated with a Green (λ G =530 nm) and a Red (λ R =627 nm) light emitting diode (LEDs). These diodes are sequentially triggered and N sums of 10 images are acquired for each illumination wavelength. Building a single sum requires about 1s.

A particular attention is paid to ensure a non-symmetric, wavelength-independent, light distribution over the pupil of the objective lens. A dichroic mirror is used to similarly illuminate a diaphragm with the two light sources. The image of the diaphragm through a collection lens is placed in the back focal plane of the objective lens (see Fig. 2).

The light distribution is made non-symmetric in the back focal plane of the objective lens by making the image of the diaphragm smaller than the pupil and by shifting this image with respect to the objective lens axis. The image of the diaphragm is thus denoted D in Fig. 3. Assuming that the observed sample is a perfectly flat mirror orthogonal to the objective lens axis, the collected light is distributed over the area C in the back focal plane. Slightly tilting the sample then shifts this area in the back focal plane. If one chooses to move the area C to the area C t (for instance), some of the reflected light is no longer collected by the objective lens (i.e., is cut out by the pupil) and the collected intensity strongly depends on the sample's surface orientation.

Reflected intensity changes may then arise either from surface reflectivity or from collection efficiency changes. The former is usually wavelength-dependent (as with electro-or thermo-reflectance [START_REF] Tessier | Measuring and predicting the thermoreflectance sensitivity as a function of wavelength on encapsulated materials[END_REF] ) while the latter is (at first glance) wavelength-independent since it primarily depends on the surface orientation. This property is thus exploited in the following to provide both local surface modification and kinematic field measurements. The corresponding formalism previously introduced 12 is first recalled. The intensity I m collected by the pixel P at wavelength λ reads

I m (P, λ, t n ) = I i (P, λ)R(P, λ, t n )ε (θ(P ), t n ) ( 1
)
where I i (P, λ) is the intensity impinging on the surface element conjugated with pixel P , at wavelength λ. R(P, λ, t n ) is the reflectivity of the surface, i.e., the ratio of the reflected and incident intensities, depending on a local driving parameter X(P, t n ) (such as electrical charge density, temperature, etc.) at the time step t n . ε (θ(P ), t n ) is the collection rate related to the local surface orientation θ(P ). The influence of X(P, t n ) on the reflectivity is assumed to be small, so that R is linearized

R(P, λ, t n ) = R 0 (P, λ)[1 + r(λ)X(P, t n )] (2)
with r(λ) = 1 R 0

∂R

∂X the relative reflectance sensitivity to the driving parameter X. The collection rate ε (θ(P ), t n ) is linearized in the case of small surface rotations around the initial orientation θ 0 (P )

ε(θ(P ), t n ) = ε 0 (θ 0 (P )) [1 + d θ 0 (P )(θ(P, t n ) -θ 0 (P ))] (3)
with d θ 0 (P ) = 1 ε 0 ∂ε ∂θ the local slope sensitivity. Finally, the measured intensity linearly depends on the wavelength-independent and -dependent relative intensity changes R wi (P, t) and R wd (λ, P, t), respectively :

I m (P, λ, t n ) = I a (P, λ) [1 + R wd (λ, P, t n ) + R wi (P, t n )] (4)
with

I a (P, λ) = I i (P, λ)R 0 (P, λ)ε 0 (θ 0 (P )) (5) R wd (λ, P, t n ) = r(λ)X(P, t n ) (6) R wi (P, t n ) = d θ 0 (P )(θ(P, t n ) -θ 0 (P )) (7)
For a given (P, t n ), Eq.( 4) taken for both analyzed wavelengths are recast as the linear system :

  I m (λ R ) -I a (λ R ) I m (λ G ) -I a (λ G )   =   I a (λ R ) I a (λ R ) k I a (λ G ) I a (λ G )     R wd (λ R ) R wi   (8)
where the ratio k = r(λ G )/r(λ R ) is assumed to be different from 1 in order to ensure the system is invertible.

III. IDENTIFICATION OF THE RELATIVE REFLECTANCE SENSITIVITIES RATIO

The knowledge of the k value corresponding to the considered physical or chemical phenomena is therefore crucial to a robust decoupling procedure. A method to attain it without any additional experimental effort and without any additional assumption is described in the following. The results illustrating the full procedure are obtained by sweeping the electrode potential U during three cycles between -0.2 and 0.2 V at 4mV.s -1 while recording the electrode charge Q(t n ) = X(P, t n )dS (dS is the elementary surface). N = 44 sums of 10 images are sequentially acquired for each wavelength during each cycle, and 10 additional sums of 10 images have been acquired before sweeping the potential. Stacks for which (at least) one image is found missing are removed, so that 103 images stacks are kept.

The recorded voltammogram is displayed on Fig. 4 and indicates that the electrochemical system stays in the double-layer regime, meaning no chemical reaction (besides adsorption) occurs at the electrode's surface. After the initial cycle, all subsequent voltage cycles lead to identical responses.

A. Identification of the k-ratio

In order to lower the experimental complexity, the proposed procedure makes use of the acquired images, only keeping the red and green intensities collected by pixels corresponding to the part of the electrode on top of the rigid substrate (Ω s ). This is assumed not to deform during the experiment, so that the recorded reflectivity changes arise only from chemical modifications of the surface. Eq.4 is then recast assuming that

R wi (P, t n ) = 0 ∀P ∈ Ω s , ∀t n (9)
so that, composing Eq. ( 4) by ln and assuming small reflectivity changes

ln (I m (P, λ, t n )) -ln (I a (P, λ)) ≃ r(λ)X(P, t n ) ∀P ∈ Ω s , ∀t n (10)
The equations for the different wavelengths are then easily combined and averaged to cancel the contribution of I a (P, λ) :

L(P ∈ Ω s , λ R ) -k -1 L(P ∈ Ω s , λ G ) = 0 (11)
where the vector L is defined, for pixels P ∈ Ω s , by

L(P, λ) = ln (I m (P, λ, t n )) -ln (I m (P, λ, t n )) (12)
and • stands for the time average. L(P, λ) therefore features one entry per pixel and per time step. As a consequence of Eq. ( 11), the matrix

M k M k = [L(P ∈ Ω s , λ R ), L(P ∈ Ω s , λ G )] (13)
should be singular. A singular value decomposition [START_REF] Demmel | Jacobi's method is more accurate than QR[END_REF] (SVD) of M t p M p is thus performed to yield

M t k M k = VEV t (14)
V is a 2 × 2 orthogonal matrix and E is a 2 × 2 diagonal matrix, with singular values decreasing with the index. The ratio k is therefore obtained as

k = - V 12 V 22 (15)
The identification quality can therefore be assessed by considering the ratio

S k = E 22 /E 11
which should vanish. The identification quality may be further assessed by rewriting Eq.( 4)

I m (P, λ R , t n ) = I a (P, λ R ) (1 + r(P, λ R )Q(t n )l(P )) ∀P ∈ Ω s , ∀t n I m (P, λ G , t n ) = I a (P, λ G ) (1 + kr(P, λ R )Q(t n )l(P )) ∀P ∈ Ω s , ∀t n (16)
with

X(P, t n ) = Q(t n )l(P ) (17)
l(P ) is a localization function relating the local charge density to the overall electrode charge. Setting k to the above-determined value, the system ( 16) may be solved in a leastsquare sense for [I a (P, λ R ), I a (P, λ G ), r(P, λ R )l(P )I a (P, λ R ), r(P, λ R )l(P )I a (P, λ G )] t at each location. One defines

η 2 r (P, t n ) = [I m (P, λ R , t n ) -I a (P, λ R ) (1 + r(P, λ R )Q(t n )l(P ))] 2 2σ 2 R + [I m (P, λ G , t n ) -I a (P, λ G ) (1 + kr(P, λ R )Q(t n )l(P ))] 2 2σ 2 G ( 18
)
where σ 2 R and σ 2 G are the estimated variances on the red and green intensity measurements, respectively. η 2 r (P, t n ) is insensitive to non-uniform charge distribution so that it is possible to further assess the identification quality by plotting the map η 2 r (P, t n ) , which can then be used to locate areas featuring a modified surface chemical activity (i.e., a different k). It should be emphasized that this identification procedure does not make use of any assumption regarding the charge distribution over Ω s (i.e., the localization function l(P )).

B. Example

This procedure is applied to the experimental results obtained during the experiment presented above. It should be highlighted that this calibration does not require any additional experimental effort, since it is based on the images acquired to monitor the coupled phenomenon under scrutiny. The only requirement is to have optical access to a fixed and rigid part undergoing the same chemical process. In the present case, Ω s covers 19 × 17 pixels. The procedure yields k = 1.18 and a ratio S k = 0.5, thereby proving the acceptable identification quality. This is also visually confirmed by plotting L(P ∈ Ω s , λ R ) as a function of L(P ∈ Ω s , λ G ), which is found to exhibit a linear correlation (see Fig. 5). The identification quality may be further assessed by mapping the local residual η 2 r (P, t n ) over Ω s (see Fig. 6). The maximum value corresponds to almost 3 times the variance on the measured intensities and thus remains moderate, thereby proving the considered assumptions are sounded.

IV. IDENTIFICATION OF THE TILT SENSITIVITY

The use of the decoupling procedure described in section II B together with the k value obtained from the identification procedure detailed in section III A allows one to compute a wavelength-dependent (R wd (λ R , P, t)) as well as a wavelength-independent (R wi (P, t)) reflectivity change fields. The former is assumed to describe the spatial distribution of the chemical modifications the surface undergoes, and the latter is assumed to reflect the mechanical deformation of the surface, as detailed in Sect.II B. A calibration procedure exploiting the bimorph effect of the structures under scrutiny is then presented in order to quantitatively 10 relate the obtained wavelength-independent reflectivity change field to the actual surface rotation field, i.e., to experimentally obtain the tilt sensitivity field d θ 0 (P ).

A. Thermal calibration procedure

The proposed calibration procedure makes use of the fact the considered structures are silicon nitride / metal beams and therefore display some bimorph effect. Controlling the temperature is thus a convenient way to prescribe a well known deformation, so that measuring the induced wavelength-independent reflectivity changes allows one to retrieve the map d θ 0 (P ) defined in Eq. ( 5). For the present case, the temperature T is varied from 24 • C to 29 • C by 0.5 • C steps, and green and red images are acquired for each step. N T = 15 images are available for each wavelength.

This temperature changes are assumed to be small enough to neglect thermoreflectance effects (∆R/R in the 10 -3 range for a 5 • C temperature change 8 ), and the temperature-induced mechanical effect is assumed to be linear with respect to the temperature change :

I m (P, λ, T n ) = I a (P, λ) [1 + d θ 0 (P )(θ(P, T n ) -θ 0 (P ))] = I a (P, λ) [1 + d θ 0 (P )βf (P )(T n -T 0 )] = I a (P, λ) 1 + ∂R wi (P, T 0 ) ∂T (T n -T 0 ) (19)
θ 0 (P ) is chosen to represent the mechanical state of the structure in any user-defined reference (temperature) state T 0 . β is a scaling factor and f (P ) is the deformation mode. As an example, a cantilever beam along the x direction (x = 0 at the base) subjected to bimetallic effect theoretically yields f (P ) = x ∀x 0 17 . Denoting M T the matrix

M T = [1, T -T 0 ] (20)
The set of Eqs 19 for all temperatures is recast

M T   I a (P, λ) I a (P, λ) ∂R wi (P,T 0 ) ∂T   = I m (P, λ) (21)
The system (21) features N T equations for two unknowns I a (P, λ) and I a (P, λ) ∂R wi (P,T 0 ) ∂T and is thus solved in a least-square sense. The presumably wavelength-independent thermal sensitivity field ∂R wi (P,T 0 ) ∂T is then retrieved by forming the ratio of the two unknowns. Fig. 7 thus displays the ∂R wi (P,T 0 ) ∂T field, averaged over the cantilever's width, as computed from green and red images, independently.

FIG. 7. ∂R wi (P,T 0 ) ∂T field, averaged over the cantilever's width, as computed from green (green line)

and red images (red line), independently.

The obtained thermal sensitivity fields share the same dependence to the position along the cantilever : the sensitivity is constant (and comparatively small) over Ω s and varies linearly along the cantilever. Additional experiments (not reported herein) allow to attribute the non-zero sensitivity over Ω s to a slight out-of-plane chip displacement induced by the temperature change of the sample holder. A slight discrepancy is also observed between the thermal sensitivities computed using green or red images, apparently violating the assumption that a purely mechanical effect induces a wavelength-independent reflectivity change. This is confirmed by plotting the so-computed red thermal sensitivity as a function of the green thermal sensitivity (Fig. 8). These are linearly related but not equal, the equality relation being plotted as a plain line in Fig. 8. Sticking to the formalism defined in Sect.

II B, one may define k

i k i = ∂R wi (P,T 0 ) ∂T G ∂R wi (P,T 0 ) ∂T R (22)
Using the SVD-based estimator described in Sect.III A, one obtains k i = 0.88 (without correcting for the non-zero sensitivity on the substrate) which has been found to be a very reproducible value for a given objective lens, so that the slight deviation from k i = 1 is attributed to chromatic aberrations. It should be highlighted that this unexpected behavior does not significantly challenge the approach presented in Sect.II B : Eq. ( 8) only has to be changed to

  I m (λ R ) -I a (λ R ) I m (λ G ) -I a (λ G )   =   I a (λ R ) I a (λ R ) k I a (λ G ) k i I a (λ G )     R wd (λ R ) R wi (λ R )   (23)
and the system is invertible if k = k i . The mechanically-induced reflectivity change and the local slope sensitivity are thus respectively denoted R wi (λ, P, t n ) and d θ 0 (λ, P ) for the sake of generality in the following.

Inverting Eq.( 7) to retrieve the surface rotation from the mechanically-induced reflectivity change requires the knowledge of d -1 θ 0 (λ, P ). It should however be highlighted that as one moves from the cantilever tip to the base, one experiences a mechanical effect of decreasing amplitude so that the reflectivity change vanishes close by the cantilever base, as seen from Fig. 7 (that is f (P ) goes to 0 when moving to the base). As a consequence, determining

d -1
θ 0 (λ, P ) from ∂R wi (λ,P,T 0 )

∂T is expected to fail in a purely local approach (see Eq. 19).

The field d -1 θ 0 (λ, P ) is thus expanded onto a user-defined basis {h n d (P )} (which is not particularized here for the sake of generality):

d -1 θ 0 (λ, P ) = N d n d =0 g n d (λ)h n d (P ) (24)
The coefficients {g n d (λ)} may thus be retrieved by solving (see Eq.( 19)) :

β -1 N d n d =0 g n d (λ)h n d (P ) ∂R wi (λ, P, T 0 ) ∂T = f (P ) (25)
Keeping {g n d (λ R )} as unknowns, one has to simultaneously solve

β -1 N d n d =0 g n d (λ R )h n d (P ) ∂R wi (λ R , P, T 0 ) ∂T = f (P ) β -1 N d n d =0 k -1 i g n d (λ R )h n d (P ) ∂R wi (λ G , P, T 0 ) ∂T = f (P ) ( 26
)
Any considered pixel provides such equalities, so that this equations set is usually overdetermined, since it features only N d + 1 unknowns. The system (26) is thus solved in a least-square sense to provide the optimal values for the unknowns β -1 g n θ (λ R ). Assuming the geometry of the cantilever is well known, β may be for instance calculated using elasticity theory [START_REF] Timoshenko | Analysis of bi-metal thermostats[END_REF] to provide the coefficients g n θ (λ R ) and thus the optimal field d -1

θ 0 (λ R , P ) = N d n d =0 g n θ (λ R )h n d (P )
. β may also be experimentally approached, as detailed in next Section.

B. Interferometric reference measurements

The margins on the dimensions of micrometer-sized devices may however be relatively large, so that this section is intended to provide an experimental procedure to retrieve β for the cantilever under scrutiny. The coefficient β is obtained prior to the above-described calibration procedure from similar experiments conducted in air with a commercially available Mirau interferometric objective lens (Nikon IC Epi Plan DI 10x). An image acquired at a given temperature T results from the phase map φ(P, λ, T ) and formally reads

I r (P, λ, T ) = I 0r (P, λ) + A(P, λ) cos [φ(P, λ, T )] ( 27
)
which is recast to highlight the role of the thermally-induced cantilever bending 17 :

I r (P, λ, T n ) = I 0r (P, λ) + A(P, λ) cos φ(P, λ, T 0 ) -φ(P ∈ Ω s , λ, T 0 ) + 2πnβ(T n -T 0 )x 2 Γλ + φ(P ∈ Ω s , λ, T n ) (28)
n is the refractive index of the surrounding medium. I 0r (P, λ) and A(P, λ) are the mean intensity and contrast fields, respectively. The Mirau objective is illuminated using the same light distribution as in Sect.II B and Γ is the correction factor introduced to account for the resulting incidence angle. Denoting

∆ s (λ, T n ) = φ(P ∈ Ω s , λ, T n ) -φ(P ∈ Ω s , λ, T 0 ) (29)
Eq.(28) reads

I r (P, λ, T n ) = I 0r (P, λ) + cos (∆ s (λ, T n )) cos 2πnβ(T n -T 0 )x 2 Γλ -sin (∆ s (λ, T n )) sin 2πnβ(T n -T 0 )x 2 Γλ A(P, λ) cos φ(P, λ, T 0 ) + -sin (∆ s (λ, T n )) cos 2πnβ(T n -T 0 )x 2 Γλ -cos (∆ s (λ, T n )) sin 2πnβ(T n -T 0 )x 2 Γλ A(P, λ) sin φ(P, λ, T 0 ) (30)
= I 0r (P, λ) +

A(P, λ) [f c (∆ s (λ, T n ), β, T n ) cos φ(P, λ, T 0 ) + f s (∆ s (λ, T n ), β, T n ) sin φ(P, λ, T 0 )] (31)
Again, for a given wavelength, the equations set (30) features N T times the number of pixels equations, for 3 unknowns per pixel (namely, I 0r , A cos φ(P, λ, T 0 ) and A sin φ(P, λ, T 0 )), in addition to the sought parameter β and to N T -1 phase shifts ∆ s (λ, T n ). The system is thus overdetermined, assuming N T is large enough. It is solved in a least-square sense, taking advantage of the structure of Eq.(30). The objective function reads η 2 (I 0r , A cos φ(P, λ, T 0 ), A sin φ(P, λ, T 0 ), β, ∆ s ) = P,Tn {I r (P, λ, T n ) -I 0r (P, λ)

-A(P, λ) [f c (∆ s (λ, T n ), β, T n ) cos φ(P, λ, T 0 ) +f s (∆ s (λ, T n ), β, T n ) sin φ(P, λ, T 0 )]} 2 (32)
The minimization problem is decomposed : min

I 0r ,A,φ(P,λ,T 0 ),β,∆s η 2 (I 0r , A, φ(P, λ, T 0 ), β, ∆ s ) = min β,∆s η 2 Ĩ 0r , Ã, φ(P, λ, T 0 ), β, ∆ s (33) = min β η 2 Ĩ 0r , Ã, φ(P, λ, T 0 ), β, ∆s (34)
β and ∆ s (λ, T n ) are first initialized, so that the optimal fields Ĩ0r , à cos φ(P, λ, T 0 ) and à sin φ(P, λ, T 0 ) are easily obtained for any chosen β and ∆ s (λ, T n ) (linear minimization step). The objective function η 2 Ĩ 0r , Ã, φ(P, λ, T 0 ), β, ∆ s thus only depends on β and ∆ s (λ, T n ), and has to be minimized with respect to these two variables sets (Eq.( 33)). This is achieved by setting the value for β and subsequently minimizing the objective function with respect to ∆ s (λ, T n ), using a conjugate-gradient algorithm [START_REF] Press | Flannery Numerical Recipes in C++[END_REF] . The resulting objective function thus only depends on β, and is minimized using the golden-search method [START_REF] Press | Flannery Numerical Recipes in C++[END_REF] , finally yielding an experimental β value. Applying this procedure to the data recorded for the considered cantilever (not shown herein) yields β = -1.49 × 10 -5 µm -1 K -1 , so that the d -1 θ 0 (λ R , P ) field may be obtained. Fig. 9 displays the obtained inverse tilt sensitivity field d -1 θ 0 (λ R , P ) along the cantilever when using a polynomial basis up to degree 2 for {h n d (x)}. This sensitivity, which is found to be rather uniform along the cantilever, is necessary to the conversion into mechanical fields of the measured reflectivity change field (see Eq.7).

V. CONVERSION FROM THE MECHANICALLY-INDUCED REFLECTIVITY CHANGE FIELD TO KINEMATIC DATA

This section thus describes a projection procedure used to convert the measured mechanicallyinduced reflectivity change field to a surface rotation field and to filter out noise. θ 0 (λ R , P ) field obtained using a polynomial basis up to degree 2 (N d = 2).

A. Conversion to kinematic data

The mechanically-induced contribution may be given a mechanical sense by making use of the d -1 θ 0 (P ) field. The cantilever displacement field w(P, t n ) at time t n may be expanded onto a user-defined functions basis {ω nw (P )}

w(P, t n ) = Nw nw=0 W nw (t n )ω nw (P ) (35)
For an Euler-Bernoulli cantilever beam, the surface rotation reads

θ(P, t n ) -θ 0 (P ) = dw(P ) dx (36) = Nw nw=0 W nw (t n ) dω nw (P ) dx (37)
so that the coefficients W nw (t n ) are to be retrieved by projecting the displacement field onto the dωn w (P ) dx basis using Eq.( 7)

Nw nw=0 W nw (t n ) dω nw (P ) dx = d -1 θ 0 (P )R wi (P, t n ) = θ(P, t n ) -θ 0 (P ) (38)
The left-hand side of Eq.( 38) is built from the user-defined basis, whereas the right-hand side is the product of the mechanically-induced field obtained in Sect.V B and of the inverse sensitivity obtained in Sect.IV A. There is such an equation for each pixel, for N w unknowns, so that Eq.( 38) is an overdetermined linear system. It is thus solved in a least-square sense to provide the projection Θ(P, t n ) = Nw nw=0 Wnw (t n ) dωn w (P ) dx of the rotation field onto the chosen basis.

It should be highlighted that this projection of the rotation field onto a functions basis is not strictly necessary (Eq.( 7) could be solved for each pixel), but is a convenient and sound way to filter out noise. In addition, N w and N d should not be chosen independently : as seen from Eq.( 38), setting N d > N w would filter out terms kept when retrieving d -1 θ 0 (P ) and thus does not make sense. The mechanically-induced R wi (λ R , P, t n ) and wavelength-dependent R wd (λ R , P, t n ) fields are obtained for the experiment described in Sect.II A and III by solving the linear system (23) for each point and time step, with the parameters identified in Sect.III, IV A. One defines the quantities Rwd (λ R , x, t n ) and Rwi (λ R , x, t n ) as the averages over the cantilever's width of the R wd (λ R , P, t n ) and R wi (λ R , P, t n ) fields, respectively, in order to allow a 2-D description. Fig. 10 displays Rwd (λ R , x, t n ), with the cantilever clamping located at x = 16 pix and the cantilever tip at x = 163 pix. Rwd (λ R , x, t n ) is mostly about 1% in magnitude, with a notice- able heterogeneity along the cantilever. It is also seen to be non-zero and electrochemicallydriven on the substrate (x < 16 pix), thereby proving its non-mechanical origin. This field actually corresponds to the reflectivity field used to monitor surface chemical affinities [START_REF] Munteanu | Mapping fluxes of radicals from the combination of electrochemical activation and optical microscopy[END_REF] , and illustrates (wavelength-dependent) electroreflectance phenomena in this particular case [START_REF] Amiot | Identification of the electroelastic coupling from full multi-physical fields measured at the micrometre scale[END_REF] . It also shows a slight charge accumulation over the cycles. Fig. 11 shows Rwi (λ R , x, t n ), which is seen to be fairly larger in magnitude (up to 10%) for the example under scrutiny, owing to the special illumination detailed in Sect.II B. It is found to have a comparatively very low value on the substrate, which does not exhibit any time dependence, thereby proving its mechanical origin. One should however notice a perturbation around x = 58 pix at t 87 , which has been found to result from a bubble in the solution traveling above the cantilever.

B. Example

The raw rotation field θ(P, t n ) -θ 0 (P ) and its projection Θ(P, t n ) onto a polynomial basis (up to degree 3 in x) are computed according to Eq. ( 38). Again, one defines the quantities θ(x, t n )-θ0 (x) and Θ(x, t n ) as the averages over the cantilever's width of the θ(P, t n )-θ 0 (P ) and Θ(P, t n ) fields, respectively. Fig. 12 shows the raw rotation field averaged over the cantilever's width θ(x, t n ) -θ0 (x). It is not computed for locations over the substrate (x < 16 pix) since the calibration procedure does not deform this surface and thus does not allow for the identification of d -1 θ 0 (λ R , P ) over the substrate (see Fig. 9). The computed θ(x, t n ) -θ0 (x) is thus arbitrarily set to 0 for x < 16 pix in Fig. 12. One should however notice the very good continuity of the rotation at FIG. 12. Raw rotation field, averaged over the cantilever's width θ(x, t n ) -θ0 (x) (rad). FIG. 13. Projected rotation field, averaged over the cantilever's width Θ(x, t n ) (rad).

x = 16 pix, reflecting the mechanical nature of R wi (λ R , P, t n ). Cycling the potential makes it bend down (away from the electrode side) when the electrode potential is increased, and deformation slightly accumulates across the potential cycles (similarly to the charge density, as seen from Fig. 10). Similar results were previously reported with a rather different measurement system 20 . Perturbations on the R wi (λ R , P, t n ) field are directly translated to perturbations on the raw rotation field, as seen from θ(x = 58 pix, t 87 ) -θ0 (x = 58 pix). This rotation field has been projected for all time steps onto a functional basis containing polynomials up to degree 3, and its average over the cantilever's width Θ(x, t n ) is displayed in Fig. 13. Most of the comments above apply to the projected rotation field. This projection procedure is particularly shown to filter out noise and perturbations, as one could notice from Θ(x = 58 pix, t 87 ), which no longer displays any perturbation. One should highlight that having the possibility to choose a functional basis allows one to test for assumptions on the rotation field or to explicitly assume (by using Lagrange multipliers) that the projected rotation field satisfies mechanical constrains (a clamping condition Θ(x = 16 pix, t n ) = 0 for instance). The projected rotation field may then be integrated if one wishes to retrieve the out-ofplane displacement field of the cantilever. The result is displayed on Fig. 14 and illustrates the robustness of proposed approach : measuring the surface rotation (instead of the displacement) avoids the use of further derivation (seen as a high-pass filter) in the mechanical exploitation of the data, and one further takes advantage of the projection step to cancel out some measurement noise.

VI. CONCLUSION

This contribution extensively describes calibration procedures and mathematical methods used to quantitatively translate the fields obtained from multiple-wavelength reflection microscopy with an original illumination device. These procedures are kept minimal by making use of the bimorph effect displayed by the gold-coated micromechanical structures extensively used for micromechanical sensors. Most of these calibration procedures may then be automated, keeping the experimental cost minimal. In addition, the decoupling procedure (besides calibration) essentially amounts to solving a linear system. The computation cost is thus also minimal so that (quasi) real-time implementation seem attainable. For systems which do not display any bimetallic effect, any convenient reproducible mechanical actuation may be used, using a similar treatment. The chemo-mechanical phenomena at the origin of the development of cantilever-based sensors may then be thoroughly and quantitatively studied using highly redundant data from a simple commercial reflection microscope featuring minor modifications.

FIG. 1 .

 1 FIG. 1. Optical view of a cantilever beam. Ω s is the portion of rigid substrate used for calibration (see Sect. III A).

FIG. 2 .

 2 FIG.2. Schematic of the multiple wavelengths reflection microscope.

FIG. 3 .

 3 FIG.3. Representative light distribution in the objective back focal plane π (see Fig.2). D is the image of the illumination diaphragm, C is the collected light distribution for a flat orthogonal mirror. C t is the collected light distribution obtained when tilting such a sample.

FIG. 4 .

 4 FIG. 4. Recorded current during the three cycles as a function of the imposed potential gap U .

FIG. 5 .

 5 FIG. 5. L(P ∈ Ω s , λ R ) as a function of L(P ∈ Ω s , λ G ) and the line corresponding to the identified k value.

5 FIG. 6 .

 56 FIG.6. Normalized residual map η 2 r (P, t n ) over the substrate P ∈ Ω s .

FIG. 8 .

 8 FIG. 8. Red thermal sensitivity as a function of the green thermal sensitivity. The solid line corresponds to k i = 1.

FIG. 9 .

 9 FIG. 9. d -1θ 0 (λ R , P ) field obtained using a polynomial basis up to degree 2 (N d = 2).

FIG. 10 .

 10 FIG. 10. Identified wavelength-dependent reflectivity change field under cyclic electrochemical actuation, averaged over the cantilever's width Rwd (λ R , x, t n ) (unitless). The cantilever position is schematically drawn on the left, and the actual potential value is recalled on top.

FIG. 11 .

 11 FIG. 11. Identified mechanically-induced reflectivity change field, averaged over the cantilever's width Rwi (λ R , x, t n ) (unitless).

FIG. 14 .

 14 FIG. 14. Displacement field, averaged over the cantilever's width, as obtained by integrating Θ(x, t n ) (in µm).

ACKNOWLEDGEMENT

This work was partially funded through the NRA grants µEcoliers ANR-08-JCJC-0088 (F.K., S.M., Y.F.) and CheMeCo ANR-11-JS09-019-01 (C.W., F.A.). This work was partly supported by the french RENATECH network and its FEMTO-ST technological facility.


I Abed

N Kacem

N Bouhaddi

M L Bouazizi

Multi-modal vibration energy harvesting approach based on nonlinear oscillator arrays under magnetic levitation

Keywords: numbers: 41.20.Gz, 05.45.-a, 46.40.Ff, 46.15.Ff, 45.10.Db

We propose a multi-modal vibration energy harvesting approach based on arrays of coupled levitated magnets. The equations of motion which include the magnetic nonlinearity and the electromagnetic damping are solved using the harmonic balance method coupled with the asymptotic numerical method. A multi-objective optimization procedure is introduced and performed using a Non-dominated Sorting Genetic Algorithm (NSGA) for the cases of small magnet arrays in order to select the optimal solutions in term of performances by bringing the eigenmodes close to each other in terms of frequencies and amplitudes. Thanks to the nonlinear coupling and the modal interactions even for only three coupled magnets, the proposed method enable harvesting the vibration energy in the operating frequency range of 4.6 -14.5 Hz, with a bandwidth of 190 % and a normalized power of 20.2 mW cm -3 g -2 .

Introduction

The natural energy sources have been frequently taken to produce inexhaustible electric energy in their local environment. The purpose is related to the reduction of power requirement and to replace a battery that has a limited lifetime, requires maintenance and cannot be used in hostile environments. Then, the mechanical structures were exposed to ambient and internal vibrations to produce energy by using adequate transducers [1] which convert the mechanical force to electric power. The accumulated energy in these cases can be stored and used. Within the emerging trend of smart systems, Vibration Energy Harvesters (VEHs) can be used in many fields such as environmental monitors, wireless sensors and medical implants [2,[START_REF] Ding | Wireless Communications[END_REF][START_REF] Wang | Environmental Science and Information Application Technology[END_REF].

In order to make VEHs usable, there are several types of transduction, where the most common transduction modes are piezoelectric [START_REF] Ottman | Power Electronics[END_REF]6,7,8,9], and electromagnetic [10,1,11,12]. The conventional linear VEHs are usually designed to be resonantly tuned to the ambient dominant frequency. They have a narrow operating bandwidth that limits their application in real-world environments where the ambient vibrations have their energy distributed over a wide spectrum of frequencies, with significant predominance of low frequency. In order to overcome this limitation, several approaches have been proposed such as resonant frequency tuning techniques, multi-modal energy harvesting and the exploitation of nonlinearities.

Resonant frequency tuning technique can be implemented in a manual way by applying a preload [13,14,[START_REF] Hu | Smart Materials and Structures[END_REF], adjusting the pre-deflection [START_REF] Loverich | Proc. SPIE[END_REF], regulating the distance between the magnets with using a spring-screw mechanism [START_REF] Challa | [END_REF] and varying the vertical relative distance to modify the attractive magnetic force [START_REF] Reissman | Proc. SPIE[END_REF] and adjusting the gravity center of the tip mass [START_REF] Reissman | Proc. SPIE[END_REF] or in a self-tuning way [START_REF] Gu | [END_REF], applying voltage to the transducer [20], switching the shunt electrical load [21] and implementing an automatic controller [22,[START_REF] Reissman | Proc. SPIE[END_REF].

The multi-modal energy harvesting technique can be implemented by exploiting multiple bending modes of a continuous beam or by exploiting a cantilever array configuration. For instance, Roundy et al. [23] proposed multiple proof masses attached on a clamped-clamped beam and the implementation of this idea was realized by Yang et al. [24] with an electromagnetic generator based on three fixed magnets along a beam vibrating at the top of three coils and operating in three modes of vibration. Sari et al. [25] studied a structure composed of different assembly generators resonating at different resonance frequencies. This structure contains 35 beams in different lengths and each beam supports a coil that moves to the top of a stationary magnet. Yang et al. [26] analyzed a VEH based on coupled flexural vibration of two elastically and electrically connected piezoelectric beams, while Eurturk et al. [27] modeled an L-shaped cantilever piezoelectric energy harvester. Tadesse et al. [28] proposed a multimodal hybrid harvester consisting of piezoelectric crystals bonded to a cantilever beam and at its tip, a permanent magnet is attached to oscillate within a stationary coil fixed to the top of the package. Arafa et al. [START_REF] Arafa | Proc. SPIE[END_REF] presented a 2-DOFs cantilever piezoelectric VEH with a dynamic magnifier consisting of a spring-mass system which is placed between the fixed end of the piezoelectric beam and the vibrating base structure, while Kim et al. [START_REF] Kim | [END_REF] suggested a VEH composed of two piezoelectric cantilevers coupled with a common proof mass and utilizing both translational and rotational degrees of freedom.

Despite the fact that multimodal techniques enable wide bandwidth energy harvesting, they require more sophisticated interface circuits than that for a singlemode harvester due to the phase difference between the output signals in array configurations [START_REF] Tang | Broadband vibration energy harvesting techniques Advances in Energy Harvesting Methods[END_REF]. Consequently, several researches have been oriented towards the study of nonlinear systems. For instance, Daqaq et al. [START_REF] Daqaq | [END_REF] reported softening frequency response characteristics in a parametrically forced piezoelectric device with structural nonlinearities. Masana and Daqaq [33] proposed a nonlinear analysis of electromechanical generator composed by piezoelectric clamped clamped beam. Mann and Sims [34] showed analytically and experimentally how magnetic levitation could be used to extend device bandwidth through a hardening response. Nevertheless, this extension is limited by dry friction dissipation phenomenon. Mahmoudi et al. [35] proposed an alternative to overcome this issue by guiding the moving magnet vertically in an elastic way by means of sandwich beams, combined with an hybrid piezoelectric and electromagnetic transductions.

In this paper, unlike classical VEHs either linear [23,25,27,[START_REF] Kim | [END_REF] and multimodal or nonlinear and mono-frequency [START_REF] Daqaq | [END_REF]34,35], to our knowledge we propose for the first time a vibration energy harvesting approach based on an array of coupled levitated magnets combining the benefits of nonlinearities and modal interactions. These benefits consist essentially in enlarging the bandwidth and increasing the harvested power. The set of coupled nonlinear equations of motion is solved using the harmonic balance method (HBM) coupled with the asymptotic numerical method (ANM) [36]. Then, in order to improve the performances of the proposed method in terms of bandwidth and harvested power, a multi-objective optimization procedure is introduced and performed using NSGA-II algorithm [START_REF] Deb | Evolutionary Computation[END_REF] for the cases of two and three coupled magnets. Comparisons have been made between uncoupled VEHs based on a single levitated magnet, two and three coupled levitated magnets and state of the art VEHs in order to emphasize the high performances of the proposed multi-modal vibration energy harvesting method.

System modeling

Design

Inspired by vibration energy harvesting based on a single levitated magnet [34], we propose an extension of magnetic levitation for a multi-degree of freedom (MDOF) VEH. Using this approach, one can take advantage of modal interactions and magnetic nonlinearities in order to obtain high performances in terms of bandwidth and harvested power. As shown in Figure 1

Magnetic forces

In order to describe the electromagnetic induction model represented in Figure 1, n + 1 reference frames have been applied. The first reference frame is fixed in space and is used to describe the motion amplitude Y 0 and excitation frequency Ω of the outer housing.

The other reference frames (designed as x 1 , x 2 ...x n-1 and x n ) describe the motion of the moving magnets which are subjected to gravitational forces ⃗ P j = M j ⃗ g and magnetic forces (Figure 2) expressed as follows [10]:

⃗ F j,j-1 = - µ 0 4π Q M j-1 Q M j (d -x j + x j-1 ) 2 ⃗ x (1) ⃗ F j,j+1 = µ 0 4π Q M j+1 Q M j (d + x j -x j+1 ) 2 ⃗ x (2)
Figure 2: Schematic of a moving magnet subject to magnetic and gravitational forces.

where Q M j is the magnetic intensity of the jth magnet, µ 0 is the magnetic permeability, d is the gap separating each two adjacent magnets, x 0 = x n+1 = Y is the excitation amplitude and g = 9.81m/s 2 is the gravitational constant. A change of variable v j = x j -Y yields the resulting restoring force applied on each moving magnet

M j : ⃗ F m j = ⃗ F j,j-1 + ⃗ F j,j+1 + ⃗ P j = ( µ 0 Q M j 4π ( Q M j-1 (d+v j-1 -v j ) 2 - Q M j+1 (d+v j -v j+1 ) 2 ) -M j g ) ⃗ x (3)
where the magnetic intensity of each magnet M j is written as Q M j = S j H cj , with S j and H cj are respectively the cross section and coercive force of the magnet M j .

At static equilibrium, the resulting restoring force verifies ⃗ F m (j) = ⃗ 0. Assuming that the equilibrium position of each moving magnet is v j = 0, we obtain the following equation:

Q M j = 4πM j gd 2 µ 0( Q M j-1 -Q M j+1 ) (4)
Substituting Equation (4) into Equation (3), yields:

⃗ F m j = ( M j gd 2 (QM j-1 -Q M j+1 ) ( Q M j-1 (d+v j-1 -v j ) 2 - Q M j+1 (d+v j -v j+1 ) 2 ) -M j g ) ⃗ x
(5)

Electromagnetic damping forces

When the device is subjected to an external mechanical vibration, each moving magnet M j oscillates around its equilibrium position and a current is induced in each coil as shown in Figure 3, resulting in the following electrical damping forces:

F e j = α j-1 i j-1 -α j i j ( 6
)
where the electromechanical coupling coefficient α j = N j Bl, with B the average magnetic field strength, N j is the number of coil turns and l is the coil length.

Figure 3: Equivalent electro-mechanical system with n moving magnets based on magnetic levitation system.

In order to overcome the issue of phase difference between the output currents in array configurations, an independent interface circuit is proposed for each current (Figure 3). By applying Kirchhoff's theory to the electrical circuits, we obtain the following equation:

R j i j = α j ( vj+1 -vj ) -r int i j (7)
Equations ( 6) and (7) give the following expression:

F e j = ce j-1 ( vj -vj-1 ) -ce j ( vj+1 -vj ) (8)
The electrical damping can be expressed as a function of the internal resistance of the coil, the resistance of the harvesting circuit and the electromechanical coupling coefficient α j [9].

c ej = α j 2 R j + r int (9)

Equations of motion

The application of Newton's first law to each moving magnet M j leads to the following equation of motion:

M j vj + c vj + F e j + F m j = -M j Ÿ ; j = 1, 2, ...n (10)
Expanding the non-linear magnetic forces written in Equation ( 3) in Taylor series up to the third order, a system of coupled nonlinear equations is obtained as follows:

M j vj + (c + c ej-1 + c ej ) vj -c ej vj+1 -c ej-1 vj-1 +k j,j+1 (v j -v j+1 ) + k j,j-1 (v j -v j-1 ) +α j,j+1 (v j -v j+1 ) 2 -α j,j-1 (v j -v j-1 ) 2 +λ j,j+1 (v j -v j+1 ) 3 + λ j,j-1 (v j -v j-1 ) 3 = -M j Ÿ ; j = 1, 2, ...n ( 11
)
This expansion is valid for magnet displacements below 50 % of the gap [35]. The 5), and assuming that

change of variable Q M j = Q ′ M j M j is substituted into Equation (
the magnetic intensities are equal (Q ′ M j-1 = Q ′ M j+1
), the linear and nonlinear stiffnesses can be written as:

k j,j+1 = 2g d M j M j+1 (M j-1 -M j+1 ) k j,j-1 = 2g d M j M j-1 (M j-1 -M j+1 ) α j,j+1 = 3 2d k j,j+1 α j,j-1 = 3 2d k j,j-1 λ j,j+1 = 2 d 2 k j,j+1 λ j,j-1 = 2 d 2 k j,j-1 (12)
Equation ( 11) can be written in its matrix form. Then, in order to ensure the symmetry of the rigidity matrix, we assume that k j,j+1 = k j+1,j which results in the following relation between the magnets in term of mass

M j-1 -M j+1 = M j -M j+2 ; j = 1, 2, ...n -1 (13)
Moreover, the positivity of the linear rigidity terms (k j,j+1 ≥ 0) is conditioned by the following inequality:

M j-1 -M j+1 ≥ 0; j = 1, 2, ...n (14)

Power transferred to the electrical circuit

The magnetic transduction is ensured by n -1 coils. The oscillations of the movable magnets cause magnetic field variations in the separation zones, which provides an induced current (Lenz's Law). The induced current can be expressed as a vibration velocity function vj (t) = V j Ω sin(Ωt). The instantaneous electrical power P (t) can be written as follows:

P (t) = Ω 2 2 ( n+1 ∑ j=1 ( c ej (V j-1 -V j ) 2 )
)

(1 -Cos(2Ωt)) ( 15
)
When the system reaches the steady state regime, the average power delivered to the electrical load takes the following form:

P m = Ω 2π ∫ 2π Ω 0 P (t)dt = Ω 2 2 ( n+1 ∑ j=1 ( c ej (V j-1 -V j ) 2 ) ) (16)
Assuming that r int << R j , the mean power P m is a good approximation of the harvested power. The proposed configuration has the advantage of avoiding the phase difference between the induced currents. However, it requires sophisticated circuits in order to deliver the energy separately.

Solving procedure

The solving procedure is purely computational and based on the harmonic balance method (HBM) coupled with the asymptotic numerical method (ANM) [36]. This technique gives the periodic solutions of a dynamical system when a control parameter is varied. But its limitation is that the equation to solve should be quadratic and it is not always easy to recast any system to a polynomial quadratic form. Once the system is transformed into a quadratic form, the HBM is used to obtain the associated nonlinear algebraic system which is solved using the ANM. The application of this solving technique to Equations ( 11) and ( 16) is detailed in the Appendix.

Results and discussion

Specifications of a nonlinear SDOF VEH

We start by calculating the response of a single nonlinear levitated magnet under harmonic base excitation shown in Figure 4. Although, previous calculations of this problem exist in the literature [34], they are either limited to slightly nonlinear oscillators and/or the dependence of the excitation amplitude on the drive frequency is not taken into account. We solve it here while addressing both issues and as a precursor to the many-oscillator case, treated in the following sections. For n=1, the equation of motion of a Single Degree-of-Freedom (SDOF) VEH based on magnetic levitation can be written as follows:

M 1 v1 + (c + c e1 + c e2 ) v1 + (k 12 + k 10 ) v 1 + (α 12 -α 10 ) v 1 2 + (λ 12 + λ 10 ) v 1 3 = -M 1 Ÿ (17)
The natural frequency of the associated linear oscillator is:

f r = 1 2π √ k 12 + k 10 M 1 (18)
Substituting Equation (12) into Equation ( 18), the natural frequency becomes:

f r = 1 π √ g (M 0 + M 2 ) 2d (M 0 -M 2 ) ( 19
)
where M 0 ≥ M 2 as given by the inequality (14). Remarkably, the natural frequency depends only on the masses of the two fixed magnets. This is a consequence of assuming equal magnetic intensities and surfaces for the three magnets. Obviously, this is difficult to obtain in practice and the frequency may slightly depend on the mass of the moving magnet. Introducing the frequency f 0 = 1 π √ g/2d and the dimensionless parameter 19) becomes: Notably, f r has a lower bound which is f 0 when µ = 0, i.e. M 2 = 0. It corresponds to the case of a single repulsive force provided by the fixed magnet on the bottom which is indispensable to counter gravity. Moreover, Equation ( 19) is incompatible with the configuration M 2 = M 0 for which the static equilibrium position of the moving magnet changes due to multistability and hence the resonance frequency expression will be modified.

µ = M 2 /M 0 , Equation (
f r = f 0 √ 1 + µ 1 -µ (20)
Substituting Equation (12) into Equation [START_REF] Challa | [END_REF] and dividing by M 1 , we obtain the following equation:

v1 + 2ξω r v1 + ω 2 r v 1 - 3g 2d 2 v 2 1 + 2 d 2 ω 2 r v 3 1 = κ cos(Ωt) (21)
where 21) is a Duffing equation with a quadratic nonlinearity under harmonic excitation. To analyze this equation of motion, we use perturbation techniques which are well adapted to small excitation and damping, typically valid in levitated magnet-based resonators [34]. To facilitate the perturbation approach, in this case the method of averaging [START_REF] Nayfeh | Introduction to Perturbation Techniques[END_REF] for its ease of use, a standard constrained coordinate transformation is introduced, as given by:

ω r = 2πf r , ξ = c + c e1 + c e2 2M 1 ω r and κ = Y 0 Ω 2 . Equation (
v 1 = A(t) cos[Ωt + β(t)] (22)
In addition, since near-resonant behavior is the principal operating regime of the proposed system, a detuning parameter, σ, is introduced, as given by:

Ω = ω r + ϵσ ( 23
)
where ϵ is the small nondimensional bookkeeping parameter. Separating the resulting equations and averaging them over the period in the t-domain results in the system's averaged equations, in terms of amplitude and phase, which are given by:

Ȧ = 1 2ω r [ κ sin β -2ξω 2 r A ] (24)
β = -σ + 1 2Aω r [ 3ω 2 r 2d 2 A 3 + κ cos β ] (25)
The steady-state motions occur when Ȧ = β = 0, which corresponds to the singular points of Equations ( 24) and (25). Thus, the frequency response equation can be written in its parametric form with respect to the phase β as follows:

σ = 3κ 2 16d 2 ξ 2 ω 3 r sin 2 β + ξω r cot β (26) A = κ 2ξω 2 r sin β ( 27
)
The critical amplitude is the oscillation amplitude A c above which bistability occurs [START_REF] Kacem | [END_REF]40,41]. Thus, A c is the transition amplitude from the linear to the nonlinear behavior. At the critical drive, the resonance curve exhibits a point of infinite slope, called the critical point. Moreover, at the same point, the phase curve also exhibits an infinite slope at the same detuning as the resonance curve itself. Mathematically, A c is defined as the oscillation amplitude for which the equation dσ dβ = 0 has a unique solution

β c = π 3 .
Thus, the critical force is deduced:

κ c = 8 √ 2 3 4 √ 3 dξ 3/2 ω 2 r ( 28
)
The critical amplitude A c is obtained by substituting Equation ( 29) into Equation ( 27)

at the point β = π/2. A c = 4 √ 2 3 4 √ 3 d √ ξ ( 29
)
Excitation frequency Ω (rad/s) Figure 6 shows the evolution of the harvester frequency response (Equations ( 26) and ( 27)) with respect to the excitation amplitude for design parameters listed in Table 1.

These parameters have been chosen in order to obtain a SDOF VEH resonating at a frequency below 10 Hz (f r = 6.5 Hz) with a low electromagnetic damping. Beyond the critical amplitude, the bandwidth definition is reconsidered and determined with respect to the highest bifurcation point. Notably, the harvester bandwidth is significantly enlarged when the moving magnet is driven in the nonlinear regime.

Multi-objective optimization of a 2-DOFs VEH

For n=2, Equation (11) leads to the following coupled equations of motion:

                                          M 1 v1 + (c + c e1 + c e2 ) v1 -c e2 v2 +k 12 (v 1 -v 2 ) + k 10 v 1 + α 12 (v 1 -v 2 ) 3 -α 10 v 1 2 + λ 12 (v 1 -v 2 ) 3 + λ 10 v 1 3        = -M 1 Ÿ        M 2 v2 + (c + c e2 + c e3 ) v2 -c e2 v1 +k 23 v 2 + k 21 (v 2 -v 1 ) + α 23 v 2 2 -α 21 (v 2 -v 1 ) 2 + λ 23 v 2 3 + λ 21 (v 2 -v 1 ) 3        = -M 2 Ÿ (30)
In this case, we seek to take advantage of the nonlinear magnetic coupling between the oscillators in order to enhance the bandwidth and the harvested power of the considered device. To do so, two objective functions g 1 and g 2 have been defined: the first one is related to the distance between the two natural frequencies of the system f 1 and f 2 while the second one concerns the intermodal distance between the harvested powers taken at the eigenfrequencies.

g 1 = f 2 -f 1 f 1 ( 31
)
g 2 = |P 1 -P 2 | P n1 (32)
Where P 1 and P n1 are the harvested and normalized powers at the frequency f 1 (Equation ( 16) for Ω = Ω 1 = 2πf 1 ) and P 2 is the harvested power at the frequency f 2 (Equation ( 16) for Ω = Ω 2 = 2πf 2 ).

Since g 1 and g 2 depend in several design parameters, as a first step, a Sensitivity Analysis (SA) has been performed in order to select unessential parameters that can be held constant. Among several SA methods, the Morris method is a specialized randomized one-factor-at-time SA design [42] and its efficiency was improved by using a Latin hypercube sampling instead of a random sampling and radial points instead of trajectories [43]. This method was used to determine the boundary between high and low influence factors among the parameters c e1 , c e2 , c e3 ,M 1 , M 2 , M 3 and d 0 . This limit is defined as the first variation is equal or greater to d max /γ, where d max is the largest distance between two consecutive maximum factors and γ = 10.

Figure 7 shows that the first objective function g 1 depends only on the magnet masses and the parameters c e1 and c e2 have a low influence on the variation of the second objective function g 2 . Consequently, c e1 and c e2 are held constant in the optimization procedure discussed below.

In order to improve the performances of the two degrees-of-freedom nonlinear VEH, the resonance peaks must be very close to each other in terms of natural frequencies and modal harvested power amplitudes. To do so, the objective functions g 1 and g 2 are used in an optimization procedure (NSGA-II) [START_REF] Deb | Evolutionary Computation[END_REF]. The two cost functions are simultaneously minimized with respect to the five retained design parameters (M 2.

1 , M 2 , M 3 , d, c e3 ). M inimize [g 1 (x), g 2 (x)] , x = [M 1 , M 2 , M 3 , d, c e3 ] T , subject to M 3 + M 0 = M 1 + M 2 , M 0 ≥ M 2 , M 1 ≥ M 3 (33)
Table 2: Design parameters of the three selected optimal solutions in Figure 8. The frequency responses in terms of amplitude, phase and harvested power are displayed in Figure 9 for design parameters of solution 1 and an excitation amplitude The frequency responses of the three selected solutions are compared in term of normalized harvested power as shown in Figure 10. The first retained solution has the highest normalized power close to 12 mW cm -3 g -2 . The specifications of the three design solutions are given in Table 3, approving the high performances of the first optimal solution. Table 3: Performances of 2-DOFs VEHs for the three selected optimal solutions in Figure 11.

Solution 1 Solution 2 Solution 3 Range of variation [ ] g 1 (cm

Y 0 =
Solution 1 2 3 f 1 (Hz)
5.2 6.1 6.5 f 2 (Hz) 8.9 10 10.6 P n (mW cm -3 g -2 ) 10. [START_REF] Wang | Environmental Science and Information Application Technology[END_REF] 8.8 9.2 BW (%) 120 85 84

Multi-objective optimization of a 3-DOFs VEH

Similarly to the case of two moving magnets, an optimization procedure is performed using NSGA-II for n = 3. It involves four objective functions: two of them concern the distance between the eigenfrequencies while the two others represent the distance between the resonance peaks in term of harvested power.

g 1 = |P 3 -P 1 | g * 1 P n1 ( 34
)
g 2 = |P 2 -P 1 | g * 2 P n1 ( 35
)
g 3 = f 3 -f 1 g * 3 f 1 ( 36
)
g 4 = f 2 -f 1 g * 4 f 1 ( 37
)
where g * i is the extrema of the function g i . The adapted Morris method was used to eliminate parameters with a low influence on the objective functions leading to four problem is written as follows:

M inimize [g 1 (x), g 2 (x), g 3 (x), g 4 (x)] , x = [M 2 , M 3 , M 4 , d] T , subject to M 1 -M 3 = M 2 -M 4 , M 0 -M 2 = M 1 -M 2 , M 1 -M 3 = M 2 -M 4 , M 0 ≥ M 2 , M 1 ≥ M 3 , M 2 ≥ M 4 (38)
The The frequency responses of the three selected solutions are compared in term of normalized harvested power as shown in Figure 13. The first retained solution has the highest normalized power close to 25 mW cm -3 g -2 . The specifications of the three design solutions are given in Table 5, approving the high performances of the first optimal solution. Table 5: Performances of 3-DOFs VEHs for the three selected optimal solutions in Figure 11.

Solution

Comparative study

In this section, comparisons between several magnetic levitation based-VEHs will be presented in order to highlight the advantages of devices having an array of coupled nonlinear oscillators in terms of bandwidth, harvested power and operating frequency range.

Nonlinear coupling benefits

The performances of two SDOF VEHs denoted V EH 1 1 and V EH 1 2 are compared to those of a 2-DOFs VEH denoted V EH 2 . The design parameters of the 2-DOFs VEH correspond to those of the first optimal solution listed in Table 2. A significant comparison is obtained if the design parameters of V EH 1 1 and

V EH 2 1 are chosen in such a way that their natural frequencies match perfectly those of V EH 2 . The latter has its first natural frequency equal to 5.2/, Hz which is higher than the lower bound limit of the natural frequency of a SDOF VEH given by f 0 = 1 π √ g/2d

(see section 3.1) if we maintain the same gap value. Consequently, the gap between 6 and the excitation amplitude is Y 0 = 3 mm.

the magnets is chosen d = 24 mm for V EH 1 1 and V EH 1 2 while the values of the other design parameters are taken from those of V EH 2 except the top magnet masses which are adjusted to ensure a perfect matching between the natural frequencies (Table 6). The frequency responses of V EH 1 1 and V EH 1 2 in terms of amplitude, phase and harvested power are displayed in Figure 14 for design parameters listed in Table 6 and an excitation amplitude Y 0 = 3 mm. It is shown that the dynamic behavior of V EH 1 1 is linear while V EH 1 2 is characterized by a hardening frequency response below the critical amplitude. The maximum of phase between the two VEHs exceeds π/2 and the overall normalized harvested power of the two SDOF VEHs is about 6.4 mW cm -3 g -2 which lower than the one of the V EH 2 as shown in Figure 15. The performances of V EH 2 are compared to the overall performances of V EH 1 1 and

V EH 2 1 in Table 7, which emphasize the advantages of the nonlinear coupling between the moving magnets. Indeed, the magnetic coupling induces nonlinearities leading to an enhancement of the frequency bandwidth via a nonlinear spring hardening effect.

Moreover, unlike the case of two SDOF harvesters, the magnetic parameters can be tuned to modify the coupling terms in such a way that the natural frequencies of the 2-DOFs VEH become closer and consequently the modal interactions increase leading to high vibration amplitudes, i.e. the harvested power increases.

Table 7: A comparison of performances between V EH 1 1 + V EH 1 2 and V EH 2 . Parameter Frequency (Hz) P n (mW cm -3 g -2 ) BW (%) V EH 1 1 + V EH 1 2 5.3-12.6 6.4 115 V EH 2 5.1-12.1 10.4 120

Modal interaction benefits

The frequency responses of V EH 2 and a 3-DOFs VEH denoted V EH 3 are compared in term of normalized harvested power as shown in Figure 16. The design parameters of the 3-DOFs VEH correspond to those of the first optimal solution listed in Table 4. The 3-DOFs VEH has the highest normalized power close to 25 mW cm -3 g -2 . The performances of both devices are given in Table 8, approving the high performances of V EH 3 compared to V EH 2 and giving credits to the modal interactions and their benefits in such systems. The remarkable specifications of the proposed devices based on arrays of levitated magnets are highlighted in Table 9 with respect to existing magnetic harvesters.

Conclusion

In this paper, we propose a multi-modal vibration energy harvesting method based on arrays of coupled levitated magnets. The equations of motion have been derived taking into account the magnetic nonlinearity and electromagnetic damping. They have been solved using the harmonic balance method coupled with the asymptotic numerical method while considering frequency dependent excitation amplitude. The case of a SDOF VEH has been analytically investigated and close form expression of the critical amplitude was extracted. Multi-objective optimization procedures are introduced and performed using NSGA-II algorithm for the cases of two and three moving magnets in order to select the ideal solutions in term of performances. The advantages of the nonlinear coupling have been pointed out by comparing the joint performances of two uncoupled VEHs and those of a 2-DOFs VEH. The latter enables harvesting the vibration energy in the frequency range 5.1 -12 Hz with a bandwidth of 120 % and a normalized harvested power of 10.4 mW cm -3 g -2 . Moreover, a 3-DOFs VEH proves to be even more performant with an operating frequency range of 4.6 -14.5 Hz, a bandwidth of 190 % and a normalized harvested power of 20.2 mW cm -3 g -2 which highlights the benefits of modal interactions within the proposed harvesting energy method.

With these promising results, the work in progress concerns the modeling of multimodal interactions for large oscillator arrays. Indeed, when the number of coupled magnets increases, the number of excited modes increases over a wide frequency range leading to ultra-wide bandwidth devices and a high harvested power. Finally, the nonlinearity can be functionalized in such periodic structures in order to generate particular collective dynamics [48] with a large number of stable multi-mode solution branches suitable for energy harvesting.

Finally, the experimental validation of the proposed multi-modal energy harvesting method is an important task that we will attempt to accomplish in the near future.

Even though the implementation of an experimental protocol is expensive and takes time, it remains one of our principal objectives.

                         vj uj 0 0 0 0 0 0                          =                          0 0 0 0 0 Y 0 cos(Ωt) 0 0                          +                          u j U l j -w j -w j,j+1 -r 1 -r 2 -r 3 P                          +                          0 0 0 0 Ω 0 0 0                          +                          0 U q j v j 2 v j v j+1 0 0 r 1 2 P q                         
m( Ż) = c(t, Ω) + l(Z) + Ωl c1 + q(Z, Z) (41) with :

U l j =            - c+c ej-1 +c ej M j u j + c ej M j u j+1 + c ej-1 M j u j-1 - k j,j+1 M j (v j -v j+1 ) - k j,j-1 M j (v j -v j-1 ) - α j,j+1
M j (w j + w j+1 -2w j,j+1 ) + α j,j-1

M j (w j + w j-1 -2w j,j-1 )

           U q j =    - λ j,j+1
M j (v j w j -v j+1 w j+1 + 3v j w j+1 -3w j v j+1 ) -λ j,j-1 M j (v j w j -v j-1 w j-1 + 3v j w j-1 -3w j v j-1 ) + r 3 r 2    P q = -r 3 and l(Z) are linear vectors and q(Z,Z) is a quadratic vector.

The harmonic balance method is applied to Equations (41) to decompose the solution Z(t) into a truncated Fourier series.

Z(t) = Z 0 + H ∑ k=1 Z c,k cos(kΩt) + H ∑ k=1 Z s,k sin(kΩt) (43)
The column vector U , with size (2H +1)N eq, where N eq is the number of equations in Equation (41), collects the components of the Fourier series as: ] T (44) Substituting Equation (44) into Equation (41), assembling the terms of the same harmonic and neglecting the highest order, we obtain the following system:

U = [ Z T 0 + Z T c,
ΩM (U ) = C + L(U ) + Q(U, U ) (45)
where M (.), C, L(.) and Q(., .) are operators which depend on m(.), c, l(.) and q(., .) as it is defined by Cochelin et al. [36]. Once the algebraic system is obtained, we solve it using a continuation technique. In order to apply the ANM, equation ( 45) is transformed into:

R (U, ω) = C + L (U ) + Q (U, U ) -ωM (U ) = 0 (46)

  , the considered device is composed of n + 2 magnets M i where i ∈ [0, n + 1]. M 0 and M n+1 are fixed with respect to a Teflon tube inside which M 1 , M 2 ...M n-1 and M n are subjected to magnetic levitation forces. All magnets are placed vertically in such a way that all opposed surfaces have the same pole and wirewound copper coils are wrapped horizontally around the separation distance between each two adjacent magnets.

Figure 1 :

 1 Figure 1: A schematic diagram of the MDOF VEH based on magnetic levitation.

6

 6

Figure 4 :

 4 Figure 4: A schematic diagram of a SDOF VEH based on magnetic levitation.

Figure 5 :

 5 Figure 5: Variation of the natural frequency of a VEH based on a single moving magnet with respect to the gap d and the mass ratio µ.

Figure 5

 5 Figure 5 displays the variation of the natural frequency of a single moving magnet

Figure 7 :

 7 Figure 7: Application of the adapted Morris method to the objective functions g 1 and g 2 .

Figure 8 :Figure 9 :

 89 Figure 8: Preferred solutions of the multi-objective optimization problem written in Equation (33).

Figure 10 :

 10 Figure 10: Forced frequency responses of the three optimal solutions for a 2-DOFs VEH in term of normalized harvested power.

retained design parameters which are M 2 , M 3 ,Figure 11 :

 2311 Figure 11: Preferred solutions of the multi-objective optimization problem written in Equation (38).

  Figure 12 for design parameters of solution 1 and an excitation amplitude Y 0 = 2.3 mm. It is shown that the two moving magnets M 2 and M 3 vibrate either in phase or with a low phase shift in the frequency range [0 20]Hz. Remarkably, the displacement of M 1 has a low phase shift with respect to M 2 and M 3 displacements, up to an excitation frequency of 6 Hz and reaches a maximum of 2π/5 for an excitation frequency equal to

Figure 12 :Figure 13 :

 1213 Figure 12: Forced frequency responses of a 3-DOFs VEH in terms of (a) amplitude and phase (b) harvested power. The design parameters correspond to those of solution 1 in Table4and the excitation amplitude is Y 0 = 2.3 mm.

Figure 14 :

 14 Figure 14: Forced frequency responses of two SDOF VEHs in terms of (a) amplitude and phase (b) overall harvested power. The design parameters are listed in Table6and the excitation amplitude is Y 0 = 3 mm.

Figure 15 :

 15 Figure15: Forced frequency responses of V EH1 1 + V EH 1 2 and V EH 2 in term of normalized harvested power.

3 Figure 16 :

 316 Figure 16: Forced frequency responses of V EH 2 and V EH 3 in term of normalized harvested power.

2 n- 1 1

 211 ∑ j=1 (c ej (w j + w j+1 -2w j,j+1 )) ...v n u 1 ...u n w 1 ...w n w 12 ...w n-1,n r 1 r 2 r 3 P ] Tis the vector of unknowns with size N eq = 4(n -2) + 3, c(t, Ω) and l c1 are constant vectors, m( Ż)

Table 1 :

 1 Forced frequency responses of the SDOF VEH with design parameters listed in Table1for three configurations: linear (Y 0 = 0.7 mm), critical (Y 0 = 2 mm), nonlinear (Y 0 = 4.8 mm). Solid lines denote stable periodic solutions and dashed lines represent unstable periodic solutions. Design parameters of the considered SDOF VEH.

	Harmonics amplitude (m)	0.005 0.01 0.015	3dB 3dB	BW L	A c	BW NL	B 2	B 1	Nonlinear Critical Linear
		20 0	30	40	50		60		70	80
	Figure 6:

Table 4 :

 4 Design parameters of the three selected optimal solutions in Figure11.

		Solution 1 Solution 2 Solution 3 Range of variation
					[ ]
	g 1	0.27	0.48	0.62	[0 -1]
	g 2	0.01	1	0.05	[0 -1]
	g 3	0.86	0.85	0.85	[0 -1]
	g 4	0.63	0.92	0.92	[0 -1]
	M 2 (g)	90	110	128	[10 -250]
	M 3 (g)	10	67	63	[10 -250]
	M 4 (g)	10	43	50	10 -250]
	d 0 (mm)	15	15	15	[5 -15]
	c e1 (Ns/m) 1.02	1.02	1.02	-
	c e2 (Ns/m) 0.76	0.76	0.76	-
	c e3 (Ns/m) 1.24	1.24	1.24	-
	c e4 (Ns/m) 1.7	1.7	1.7	-
	c (Ns/m)	0.119	0.119	0.119	-

Table 6

 6

		: Design parameters of V EH 1 1 and V EH 1 2 .
	Parameter M 0	M 1	M 2	c e1	c e2	c	d
		(g)	(g)	(g)	(Ns/m)	(Ns/m)	(Ns/m)	(mm)
	V EH 1 1	80	80	10	1.19	0.13	0.119 24
	V EH 1 2	80	32	46	1.19	0.27	0.119 24

Table 8 :

 8 A comparison of performances between V EH 2 and V EH 3 .

	Parameter	Frequency (Hz) P n (mW cm -3 g -2 ) BW (%)
	V EH 2	5.1 -12.6	12.1	120
	V EH 3	4.6 -14.5	20.2	190

Table 9 :

 9 A comparison of performances between several energy harvesters based on magnetic levitation and the proposed VEHs.

	Refs	Frequency	Acceleration	Power	Normalized Power
		(Hz)	(g)	(mW )	(mW cm -3 g -2 )
	Abu Riduan et al. [10] 7-10	0.5	2.090	0.1
	Daniel et al. [44]	13	0.25	5.9	15.33
	Marin et al. [45]	47.75-51.75 0.2	25.5	0.54
	Cepnik et al. [46]	50	1	20.6	1.01
	Byung-chul et al. [47] 16	0.2	1.52	1.07
	V EH 2	5.1-12	0.22	335	10.4
	V EH 3	4.6-14.5	0.17	795	20.2

Acknowledgments

This project has been performed in cooperation with the Labex ACTION program (contract ANR-11-LABX-01-01).

Appendix

Equations ( 11) and ( 16) can be written in the following form:

Where w is a vector of unknowns, f is periodic in t and (Ω) is a real parameter.

The key point of this method lies in the quadratic recast of Equation ( 39) by introducing the following set of auxiliary variables:

Then, Equation ( 39) can be written in matrix form as follows :