Thomas Huguet
email: thomas.huguet@univ-smb.fr

Adrien Badel
email: adrien.badel@univ-smb.fr

Mickaël Lallart
email: mickael.lallart@insa-lyon.fr

Parametric analysis for optimized piezoelectric bistable vibration energy harvesters

This study focuses on vibration energy harvesting with piezoelectric bistable inertial generators in order to provide an alternative or a support to chemical batteries in the power supply of isolated wireless sensors (leading to a better autonomy of these devices). The objective of the present work is to provide guidelines to optimize this kind of generator through the investigation of the influence of different parameters (load resistance, mass, stiffness, buckling level) on its frequency response. These guidelines will be obtained exploiting a new analytical model for piezoelectric bistable harvesters which will be constructed based on a recent model introduced for electromagnetic bistable harvesters. This model includes the study of subharmonic behaviors which can be used to enhance the global bandwidth of the harvester as well as a stability robustness criterion which allows a better prediction of experimental observations. The new model will be validated by experimental data and finally exploited to provide guidelines for piezoelectric bistable harvester optimization.

I. Introduction

This study is part of the current research concerning vibration energy harvesting with the aim of designing generators able to power wireless electronic devices thanks to the mechanical energy present in their environment [START_REF] Cook-Chennault | Powering MEMS portable devices: a review of nonregenerative and regenerative power supply systems with special emphasis on piezoelectric energy harvesting systems[END_REF]. These wireless devices usually work with chemical batteries which present a major limitation: the amount of energy embedded with the wireless device is limited which leads to several operations during its lifetime to change the empty batteries. This point is moreover intensified with the phenomenon of self-discharge of primary batteries especially in high temperature environment for instance. The use of vibration energy harvesters as an alternative to conventional batteries (or as a support to overcome their self-discharge) would therefore allow to increase the autonomy of wireless electronic devices converting continuously the ambient vibration to power them during their entire lifetime.

Several kinds of harvesters have been studied in the last twenty years [START_REF] Siang | Review of vibration-based energy harvesting technology : Mechanism and architectural approach[END_REF]. Among them, inertial harvesters have the convenience of being fixable on a single point of the vibration source. Linear inertial harvesters (springmass) show interesting performance which nevertheless drastically collapse if the vibration frequency slightly moves from its resonance frequency. In order to increase the bandwidth proposed by vibration energy harvesters, part of the scientific community focused on nonlinear inertial harvesters [START_REF] Daqaq | On the role of nonlinearities in vibratory energy harvesting: A critical review and discussion[END_REF] and more particularly on a promising architecture: the bistable harvester [START_REF] Harne | A review of the recent research on vibration energy harvesting via bistable systems[END_REF]. The study of bistable designs is still a topical issue in the community of energy harvesting [START_REF] Wang | Dynamic and energetic characteristics of a bistable piezoelectric vibration energy harvester with an elastic magnifier[END_REF][START_REF] Liu | Stochastic response of bistable vibration energy harvesting system subject to filtered gaussian white noise[END_REF].

In order to optimize the performances of bistable harvesters, it is important to understand the influence of its different parameters on its frequency response as it was already proposed for linear harvesters [START_REF] Williams | Analysis Of A Microelectric Generator For Microsystems[END_REF]. The goal of this study is to develop and exploit a new analytical model for piezoelectric bistable harvesters in order to propose guidelines for their optimization. This model includes two new features: the study of subharmonic behaviors for which the mass oscillates periodically at a frequency n times lower than the excitation and the stability robustness criterion which characterizes the capability of the different high orbits (inter-well motions) to handle big disturbances without falling on a low orbit (disturbances that can be found in real conditions). Those features appear to be relevant for optimization issues as the first one (subharmonic behaviors) can be used to enhance the global bandwidth of bistable harvesters [START_REF] Huguet | Exploiting bistable oscillator subharmonics for magnified broadband vibration energy harvesting[END_REF] and the second one (stability robustness) allows better predictions of experimental measurements and observations [START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF]. The contribution of these two elements are illustrated in Figure 1 for an electromagnetic bistable harvester.

A few studies have already focused on optimization analysis for bistable harvesters but did not take into account neither the subharmonic behaviors nor the stability robustness criterion. Therefore, such results, although interesting, can be improved with the new model. Among these studies, Sneller et al. [START_REF] Sneller | Experimental investigation of a post-buckled piezoelectric beam with an attached central mass used to harvest energy[END_REF] highlighted the influence of the mass on the behavior of bistable harvesters: increasing the applied mass was found to broaden the response of the system and to lower the threshold for the onset of snap-through behaviour. Then Stanton et al. [START_REF] Stanton | Harmonic balance analysis of the bistable piezoelectric inertial generator[END_REF] focused on the mechanical and electrical damping: "there also exists an optimal impedance load and coupling for extracting maximum power for both interwell and intrawell oscillations". Finally, Panyam et al. [START_REF] Panyam | On approximating the effective bandwidth of bi-stable energy harvesters[END_REF] studied the influence of the load (time constant ratio), of the electromechanical coupling factor and of the buckling level (potential function): "while the time constant ratio has very little influence on the effective bandwidth of the harvester, increasing the electromechanical coupling and/or designing the potential function with deeper potential wells serve to shrink the effective bandwidth for a given level of excitation".

The present study is decomposed as follows. The first section proposes to establish the new analytical model for piezoelectric bistable harvesters based on a model recently proposed for electromagnetic bistable harvesters [START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF]. Then the second section introduces experimental data to validate this model. The third section exploits it to evaluate the influence of different parameters (load resistance, mass, stiffness, buckling level) on the frequency response of bistable harvesters. The conclusion section will finally summarize this study in the form of guidelines usable to optimize future piezoelectric bistable harvesters.

II. Model construction

The purpose of this section is to construct the analytical model for piezoelectric bistable harvester able to predict all its possible steady-state behaviors as a function of the excitation frequency (i.e., its frequency response). This model will be exploited later in this study to investigate the influence of the different parameters of the bistable harvester (mass, stiffness, buckling level, etc.) on its frequency response.

The bistable harvesters considered in this study is ob-tained with a buckled beam and a piezoelectric converter. The architecture is shown in Figure 2. It corresponds to a simplified schematic of the experimental prototype. This architecture includes a bistable oscillator with one degree of freedom and only one piezoelectric converter (in order to facilitate the assembly process). Contrary to the buckled beam bistable harvester introduced by Liu et al. [START_REF] Liu | Novel piezoelectric bistable oscillator architecture for wideband vibration energy harvesting[END_REF], the piezoelectric converter is here located on the side of the oscillator rather than along the moving beams. This choice has been preferred to reduce its displacement that can damage the electrical connections while its deformations remain the same. The beams on either side of the mass are doubled to remove one unwanted degree of freedom: the rotation of the mass around the axis perpendicular to the plane of the diagram. The global stiffness k shown in Figure 2 represents all the components of the bistable harvester that can store some elastic energy. Their respective stiffness are expressed in the same point to give the global stiffness k. The use of this global equivalent stiffness will lead to the generic Duffing type equations of piezoelectric bistable harvesters introduced later in Equation ( 5). Hence, the analysis which follows is valid for all the piezoelectric bistable harvesters governed by a Duffing type equation (a majority in the literature [START_REF] Harne | A review of the recent research on vibration energy harvesting via bistable systems[END_REF]) and does not depend on the particular architecture shown in Figure 2 which is given here as an example. It is noteworthy that the global stiffness k can be easily identified experimentally knowing that bistable harvesters behave like linear harvesters for small perturbations. The global stiffness can then be determined through the identification of the resonance frequency of the bistable device for small perturbations as proposed in the experimental section. The mass M represents the proof mass of the bistable harvester. The dynamic mass of the moving beams is neglected. The analytical model will focus on periodic steadystate behaviors of bistable harvester that are interesting for energy harvesting for a sinusoidal excitation: the common harmonic 1 high orbit for which the mass oscillates from one stable position to the other at the same frequency as the excitation and the subharmonic high orbits for which the mass oscillates from one stable position to the other at a frequency n times lower than the excitation for subharmonic n. The harmonic 1 high orbit indeed presents the working point with the maximum harvested power and the subharmonic high orbits allow to increase the global bandwidth of the bistable harvester as they exist on frequency ranges on which the harmonic 1 high orbit is not reachable [START_REF] Huguet | Exploiting bistable oscillator subharmonics for magnified broadband vibration energy harvesting[END_REF][START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF].

The mathematical model will also characterize these steady-state behaviors with two criteria: their stability to small disturbances and their stability robustness. The first criterion, common in the study of dynamic systems, will split the behaviors in two categories by analyzing their response when a small disturbance is added: stable if maintained, unstable if not. The second criterion is calculated for stable high orbits and defines their ability to handle bigger disturbances without falling on a low orbit (disturbances that can be found in real conditions) [START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF]. The more robust the high orbit, the more capable to handle external disturbances without falling on a low orbit and therefore the easier to maintain over time in real conditions. A threshold will be experimentally defined for this criterion above which the high orbits will be considered as robust enough to be used for energy harvesting in real conditions.

The following construction is based on a model that have been recently introduced for electromagnetic bistable harvester [START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF]. It is here adapted to piezoelectric bistable harvester adding a new degree of freedom brought by the capacitance of the piezoelectric converter. The construction is here shorten only focusing on the differences with the previous model. For more details the reader may refer to the construction of this previous model [START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF].

II.1. Formulation

The mechanical equation of the bistable harvester (shown in Figure 2) is obtained analyzing the different forces applied on the mass: the absolute inertial force, the damping force (linear dependence on x), the stiffness force and the force generated by the electromechanical conversion. The last two forces are projected on the vertical direction multiplying by x/l and multiplied by a factor two to take into account the two forces generated by the beams on both sides of the mass.

M( ẍ + A cos(ωt)) = -µ ẋ -2k 2L -2 l 2 -x 2 piezo length variation x l -2αv x l (1)
where x represents the position of the mass over time, A and ω are the amplitude and the frequency of the excitation, µ the mechanical damping coefficient and l the length of each of the four beams linked to the mass. The constant α represents the force factor of the piezoelectric converter and v its voltage. The other variables are shown in Figure 2. In this study, the buckling coefficient x 0 /L is small compared to one. Considering this assumption and l 2 = x 2 0 + L 2 , the mechanical equation becomes a Duffing-type equation:

M( ẍ + A cos(ωt)) ≈ 2k L 2 x 2 0 -x 2 x -µ ẋ -2αv x L (2)
Part of the mechanical power of the bistable harvester is then converted in electrical power by the piezoelectric converter. In this study, this electrical power is entirely dissipated in a resistance R directly connected to the piezoelectric element terminals. The power harvested by the device thus corresponds to the electrical power dissipated in this resistance. The electrical equation of the bistable harvester is obtained describing the electrical equilibrium between the piezoelectric converter (equivalent to a current generator proportional to the speed variation of the piezo length in parallel with a capacitance C 0 ) and the load resistance R:

C 0 v = α d dt 2L -2 l 2 -x 2 piezo length variation - v R (3)
Considering that the buckling coefficient x 0 /L is small compared to one and that l 2 = x 2 0 + L 2 , the electrical equation becomes:

C 0 v ≈ 2 α L x ẋ - v R (4)
The general equation set governing the behavior of the piezoelectric bistable harvester with low buckling level is therefore composed by Equation ( 2) and (4):

       ẍ + ω 2 0 2 x 2 x 2 0 -1 x + ω 0 Q ẋ + 2α ML xv = -A cos(ωt) C 0 v = 2 α L x ẋ - v R (5) with          ω 0 = x 0 L 4k M Q = x 0 L √ 4kM µ (6)
The constant ω 0 and Q are the natural angular frequency and the mechanical quality factor respectively of the equivalent linear harvester which is obtained for small oscillations of the mass of the bistable harvester around one of its stable position (x = x 0 + ∆x with ∆x << x 0 ) [START_REF] Liu | Novel piezoelectric bistable oscillator architecture for wideband vibration energy harvesting[END_REF]. The power harvested by the bistable harvester is equivalent to the power dissipated in the load resistance:

P harvested (t) = v(t) 2 R (7)
P mean harvested = 1 T T 0 v(t) 2 R dt (8)
In order to obtain the frequency response of the bistable harvester, it is now necessary to find all the steady-state behaviors which are solutions of Equation [START_REF] Wang | Dynamic and energetic characteristics of a bistable piezoelectric vibration energy harvester with an elastic magnifier[END_REF].

II.2. Harmonic balance

The focus is on behaviors useful for energy harvesting: the common harmonic 1 high orbit which presents the working point with the maximum harvested power and the subharmonic high orbits which allow to increase the global bandwidth of the bistable harvester as they exist on frequency ranges on which the harmonic 1 high orbit is not reachable [START_REF] Huguet | Exploiting bistable oscillator subharmonics for magnified broadband vibration energy harvesting[END_REF][START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF]. As a reminder, subharmonic n orbit is defined as a behavior for which the mass oscillates at frequency n times lower than the excitation frequency. Examples of these behaviors are shown in Figure 3 (time signals), Figure 4 (phase portraits) and Figure 5 (Fourier series decomposition). These examples have been obtained with a numerical integration of Equation ( 5) for different initial conditions which are detailed in appendix. The values of the bistable harvester parameters correspond to the prototype introduced in the experimental section and are summarized in Table 2.

It is interesting to note on those figures that it exists a factor two between the fundamental frequency of the mass displacement and the preponderant frequency of the voltage. This phenomenon is caused by the symmetry of the bistable harvester with respect to the axes x = 0: the piezoelectric converter deforms in the same way for positive x and for negative x. Thus, the voltage maxima (corresponding to the elongation maxima of the piezoelectric converter) are reached for all the displacement extrema of the mass (maxima and minima). The voltage minima (corresponding to the compression maxima of the piezoelectric converter) are reached each time the mass passes through x = 0.

It is also worth noting on Figure 5 that the mass displacement of the subharmonic high orbits with an odd order (1, 3 and 5) is only composed of odd harmonics of the fundamental frequency (h3, h5, h7...). For these behaviors, the voltage is composed of a fundamental frequency twice higher than the mass displacement and all its harmonics. The mass displacement of the subharmonic high orbits with an even order (2 and 4) is composed of all the harmonics of the fundamental frequency and a significant continuous component. This time, the voltage of these behaviors is composed of the same harmonics as the mass displacement (with a preponderant frequency twice higher than the fundamental frequency of the mass displacement). This difference between odd and even subharmonic orbits is directly linked to the bistability of the harvester and especially to the force applied by the spring to the mass which can be expressed as a polynomial equation involving only odd terms Ax 3 + Bx with (A, B) constants (Equation ( 5)). If f is the fundamental frequency of the mass position, the spring force only creates odd harmonics of f when the continuous component of the mass position is zero and creates all the harmonics of f otherwise. However, a behavior can only be maintained in a steady-state if the excitation can bring energy to the harvester by interacting with one of the harmonics of the mass position. For odd subharmonic behaviors, the excitation must interact with an odd harmonic of the mass position (see Figure 5). For even subharmonic behaviors, it must interact with an even harmonic of the mass position (see Figure 5). Thus, only the even subharmonic behaviors have a continuous component on the mass position which implies that their whole displacement oscillations are shifted up or down and are no longer symmetric with respect to the axes x = 0 contrary to odd subharmonic behaviors. The piezoelectric voltage is linked to a term A ẋx with A constant (Equation ( 5)). Therefore, its fundamental component is 2 f and it is only composed of even harmonics of f for odd subharmonic behaviors and of all the harmonics of f for even subharmonic behaviors.

The general form of the solution (corresponding to subharmonic n behavior) is thus approached by two Fourier series truncated to order N in order to be able to apply the classical method of harmonic balance [START_REF] Harne | A review of the recent research on vibration energy harvesting via bistable systems[END_REF]. Each of them is therefore composed of a fundamental component with an angular frequency ω/n followed by all its harmonics hω/n with h ∈ [[1, N]] and constant terms a 0 and p 0 :

           x(t) = a 0 + N ∑ h=1 a h cos h ω n t + b h sin h ω n t v(t) = p 0 + N ∑ h=1 p h cos h ω n t + q h sin h ω n t (9)
where the constant terms a 0 and p 0 corresponds respectively to the time averages of mass displacement and voltage at the resistance terminals. If a 0 is zero, the solution is a high orbit. If a 0 is close to x 0 or -x 0 , the solution is a low orbit. Since the desired response is in steady-state, the 4N

+ 2 unknown (a 0 , a h , b h , p 0 , p h , q h ) , h ∈ [[1, N]] are independent of time.
Equation ( 9) is injected in Equation ( 5) governing the behavior of bistable harvester. An approximated solution of the system can then be found thanks to the harmonic balance method [START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF]. Each approximated solution (x s , v s ) found with the harmonic balance method is defined by a unique set (a s0 , a sh , b sh , p s0 , p sh , q sh ) with h ∈ [ [1, N]]. These unique sets all correspond to possible steady-state subharmonic n behaviors of the bistable harvester for an excitation A cos(ωt).

The frequency response of the bistable harvester can therefore be obtained solving Equation [START_REF] Wang | Dynamic and energetic characteristics of a bistable piezoelectric vibration energy harvester with an elastic magnifier[END_REF] with this method for different values of n (giving the different subharmonic n behaviors including harmonic 1 for n = 1) and for different angular frequencies of the excitation ω. This method has often been applied in the literature to study the harmonic 1 behavior with only one term in the truncated Fourier series (N = 1). In our case, N has been chosen bigger than one to reduce the error on the solution sought and to reveal the subharmonic behaviors. It is determined with the following method: N is increased until the difference between the solution found and the solution for N -1 becomes negligible. This method converged to N = 6.

II.3. Stability to small disturbances

Stability to small disturbances is then evaluated for all steady-state behaviors (x s , v s ), corresponding to approximated solutions of Equation [START_REF] Wang | Dynamic and energetic characteristics of a bistable piezoelectric vibration energy harvester with an elastic magnifier[END_REF]. A small disturbance (z, w) is therefore added to the studied behavior as a pulse at the instant t = t 0 :

x(t) = x s (t) + z(t) v(t) = v s (t) + w(t) (10) with ||z(t)|| << ||x s (t)|| ||w(t)|| << ||v s (t)|| (11)
and

pulse =   z(t 0 ) ż(t 0 ) w(t 0 )   (12)
The solution (x s , v s ) is stable to small disturbances if and only if (z, w) tend towards 0. Equation ( 10) is injected in Equation ( 5) governing the behavior of the bistable harvester and linearized with respect to (z, w) to give the following differential equation governing the behavior of the disturbance:

  ż(t) z(t) ẇ(t)   = M(t)   z(t) ż(t) w(t)   (13)
with

M(t) = (14)       0 1 0 ω 2 0 2 - 3ω 2 0 2x 2 0 x s (t) 2 - 2α ML v s (t) - ω 0 Q - 2α ML x s (t) 2α C 0 L ẋs (t) 2α C 0 L x s (t) - 1 C 0 R      
The stability to small disturbances can now be evaluated with the Floquet theory and the Lyapunov exponents.

The stability of the solutions (x s , v s ) is given by the eigenvalues (β 1 , β 2 , β 3 ) of the resolvent matrix R(2T, 0), obtained after two numerical integrations of Equation ( 13) [START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF]:

R(2T, 0) =     z 1 (2T) ż1 (2T) w 1 (2T)     z 2 (2T) ż2 (2T) w 2 (2T)     z 3 (2T) ż3 (2T) w 3 (2T)     (15)
with

                                 z 1 (0) ż1 (0) w 1 (0)   =   1 0 0   numerical integration of --------------→ Equation (13) from 0 to 2T   z 1 (2T) ż1 (2T) w 1 (2T)     z 2 (0) ż2 (0) w 2 (0)   =   0 1 0   numerical intégration --------------→ Equation (13) from 0 to 2T   z 2 (2T) ż2 (2T) w 2 (2T)     z 3 (0) ż3 (0) w 3 (0)   =   0 0 1   numerical intégration --------------→ Equation (13) from 0 to 2T   z 3 (2T) ż3 (2T) w 3 (2T)   (16) and (x s , v s ) stable ⇔      |β 1 | < 1 |β 2 | < 1 |β 3 | < 1 (17)

II.4. Stability robustness

The stability robustness will now be evaluated for all the high orbits defined as stable to small disturbances. This criterion will define the ability of these high orbits to handle bigger disturbances without falling on a low orbit (disturbances that can be found in real conditions). The more robust the high orbit, the more capable to handle external disturbances without falling on a low orbit and therefore the easier to maintain over time in real conditions. A threshold will be defined for this criterion above which the high orbits will be considered as robust enough to be used for energy harvesting in real conditions.

Hence, a disturbance (z, w) is added to the stable (to small disturbances) high orbit (x s , v s ) as a speed pulse Λ at the instant t = t 0 . This pulse is not considered as negligible anymore:

x(t) = x s (t) + z(t) v(t) = v s (t) + w(t) (18)
and

pulse =   z(t 0 ) ż(t 0 ) w(t 0 )   =   0 Λ 0   (19)
The stability robustness is evaluated calculating the minimum value of the pulse Λ which leads to a destabilization of the high orbit. Thus, the high orbit remains stable (z and w tend towards 0) if Λ < Λ min and the high orbit becomes unstable (z and w tend towards infinity) if Λ > Λ min . More precisely, the stability robustness is defined by the ratio between the minimum energy needed to destabilize the high orbit and the average amount of energy provided by the excitation during one period (normative value):

Robustness = 1 2 MΛ 2 min E excitation (20)
with

E excitation = - 1 n T 0 MA cos(ωt) ẋs dt (21)
The next step is to determine the expression of Λ min . Equation ( 18) is injected in Equation ( 5) governing the behavior of the bistable harvester to give the following differential equation governing the behavior of the disturbance:

  ż(t) z(t) ẇ(t)   = M(t)   z(t) ż(t) w(t)   linear term (22) +       0 - 3ω 2 0 2x 2 0 x s (t)z(t) 2 - ω 2 0 2x 2 0 z(t) 3 - 2α ML z(t)w(t) 2α C 0 L z(t) ż(t)       non linear term
The matrix M(t) is identical to the one introduced in the small disturbance analyses in Equation ( 14). The minimal speed pulse Λ min which leads to the divergence of the disturbance (z, w) can be calculated thanks to a reductio ad absurdum applied on this differential equation. The expression of the minimal pulse Λ min which leads to destabilize the high orbit ends up to be [START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF]:

Λ min = δ 0 -Bδ 2 0 -Dδ 3 0 C (1,2) + 2α C 0 Lλ δ 0 C (1,3) (23)
with

                               δ 0 = -B + √ B 2 + 3D 3D B = 3ω 2 0 2x 2 0 max t 0 ∈[0,T] |Ψ| + 2α MLλ C (1,2) + 2α C 0 Lλ C (1,3) D = ω 2 0 4λx 2 0 C (1,2) C (i,j) = max t 0 ∈[0,T] t∈[t 0 ,t 0 +2T] e λ(t-t 0 ) R(t, t 0 ) (i,j) (24) and                                λ = - max((ln(β 1 ), ln(β 2 ), ln(β 3 )) 2T Ψ = N ∑ k=1 a sk λ 2 + ω 2 k (λ cos(ω k t 0 ) -ω k sin(ω k t 0 )) + N ∑ k=1 b sk λ 2 + ω 2 k (λ sin(ω k t 0 ) + ω k cos(ω k t 0 )) T = 2πn/ω ω k = kω/n (25)
The stability robustness can now be calculated for all the high orbits defined as stable to small disturbances knowing the minimal pulse Λ min which leads to their fall on a low orbit.

The mathematical model able to predict the frequency response of piezoelectric bistable harvesters is now complete. The next section provides experimental data to validate it. All the results (theoretical and experimental) are presented at the end of the next section in Figure 9.

III. Model validation with experimental data III.1. Prototype and identification

The experimental analysis conducted to validate the model built in the previous section has been done with the prototype shown in Figure 6. The actuator visible on the right of the top picture is not used in this study. Its stiffness is high compared to the stiffness of the piezoelectric converter. As long as this actuator is not powered, it thus does not affect the behavior of the bistable harvester.

The prototype is buckled thanks to a nut which (when it is screwed on) causes a translation of the screw linked to the actuator. The beams on both sides of the mass are doubled in order to prevent the mass from rotating.

Technically speaking, the prototype was obtained by cutting a block of APX4 steel by electrical discharge machining to make the mass, the beams and the frame appear. The piezoelectric converter (and the actuator) is then assembled. The piezoelectric converter is manufactured by Cedrat Technologies under the reference APA 120S (Amplified Piezoelectric Actuator). It consists of two elements: an elliptical stainless steel shell and a stack of inter-digitated piezoelectric ceramics placed on the long axis of this shell. The elliptical shape of the shell acts as a mechanical transformer. It amplifies the stresses applied to its small axis along its long axis. On the other hand, it reduces the deformations applied to its small axis along its long axis. This APA 120S thus optimizes the use of piezoelectric ceramics, which can support high forces but small deformations. The main dimensions of the prototype are detailed in Table 1.

This prototype corresponds to the simplified schematic introduced in Figure 2. The pivoting links are ensured by the flexibility of the thin steel beams (flexible connections). The global stiffness k shown in Figure 2 thus includes, in this particular case, (i) the stiffness of the piezoelectric converter and (ii) the stiffness of the flexible beams. The latter is here negligible as the potential energy stored in the flexible bonds is small compared to the potential energy stored in the piezoelectric converter. The stiffness of the piezoelectric converter is itself composed by the stiffness of the piezoelectric stack, the stiffness of the elliptical shell and it also includes the effect of the mechanical transformation caused by the shell which amplifies the stresses applied to its small axis along its long axis and reduces the deformations applied to its small axis along its long axis. This transformation therefore affects the global stiffness of the piezoelectric converter.

These flexible bonds have a negligible effect on the stiffness and damping of the bistable harvester and can therefore be modeled as pure pivoting links. Indeed, (i) the potential energy stored in these flexible bonds is small compared to the potential energy stored in the Piezo harvester (for low buckling configuration) and (ii) the mechanical losses in these bonds are negligible thanks to the material used (APX4 steel).

The model parameters (ω 0 , Q, C 0 , α) have been identified measuring the complex admittance Y of the piezoelectric converter (mounted on the prototype) around the angular frequency ω 0 of the bistable harvester without any excitation [START_REF] Badel | Nonlinear Conditioning Circuits for Piezoelectric Energy Harvesters[END_REF]. The different parameters are then identified adjusting the experimental admittance to the theoretical one obtained for a small voltage excitation of the piezoelectric converter inducing small oscillations of the mass around one of its stable positions: x = x 0 + ∆x with ∆x << x 0 . This identification is shown in Figure 7. The expression of the theoretical admittance being for the bistable harvester:

                 Y = I V = jC 0 ω     1 + k 2 em 1 -(ω/ω 0 ) 2 + j ω/ω 0 Q     k 2 em = 4 x 0 L 2 α 2 MC 0 ω 2 0 (26)
The other model parameters (M, L, x 0 ) are directly measured on the prototype. All the parameters and their identified values for the prototype are summarized in Table 2.

III.2. Experimental results

The frequency response obtained with the mathematical model is here compared with the frequency response obtained experimentally in order to control its relevance.

This experimental frequency response is determined as follows. The prototype is placed on a shaker ensuring a sinusoidal acceleration of constant amplitude 2.5 m.s -2 . The position of the mass and the position of the frame are measured using two laser vibrometers. The frequency ranges of harmonic 1 and subharmonic behaviors of the bistable harvester are determined with the following method: (i) The frequency of the excitation is set to a value where the behavior under study exists (based on theoretical predictions); (ii) Square voltage pulses are applied to the piezoelectric converter (30 V for 20 ms)thus generating force pulses on the mass. Each force pulse disturbs the current steady-state behavior and are likely to lead to another steady-state behavior (orbit jump), including possible subharmonics. These pulses are applied until the studied behavior is reached; (iii) Once the desired behavior is reached, the excitation frequency is slowly and smoothly increased or decreased in order to explore its entire frequency range without destabilizing it. The experimental setup is shown in Figure 8.

Figure 9 shows the theoretical and experimental fre- quency response of the piezoelectric bistable harvester. Only the theoretical behaviors defined as stable to small disturbances are drawn on this figure (the unstable behaviors are not interesting for energy harvesting as they cannot be reached experimentally). The stable behaviors are split in two different categories depending on the stability robustness criterion: robust enough or not robust enough to be maintained over time in real conditions. Each theoretical and experimental behavior are characterized by four items: the amplitude of the mass oscillation (greater than x 0 for high orbits and lower than x 0 for low orbits), the phase of the excitation displacement when the position of the mass reaches a maximum (for high orbits), the stability robustness (for high orbits) and finally the mean harvested power (P mean = v 2 rms /R). The stability robustness threshold has been chosen in such a way that the cutting frequencies of the different high orbits match with the experimental observations. The phase of the excitation displacement when the position of the mass reaches a maximum (equivalent to the phase shift for harmonic 1 behavior) is calculated considering that the mass displacement signal can be approached by its fundamental frequency. In other words, it is supposed that the instant on which the mass reaches a maximum does not change if we only take the fundamental frequency. This assumption is illustrated with few examples shown in Figure 10 (numerical study). For subharmonic n high orbit, the fundamental component of the displacement is:

x fundamental (t) = a 1 cos ω n t + b 1 sin ω n t = a 2 1 + b 2 1 cos ω n t + ψ (27) with ψ = tan -1 -b 1 a 1 + π (if a 1 < 0) (28)
The instant on which the mass position reaches a maximum is called t 1 . Hence, the phase of the mass position at t 1 is equal to zero:

ω n t 1 + ψ = 0 mod[2π] ( 29
) ⇒ t 1 = -ψ n ω mod 2π n ω ( 30
)
The excitation displacement can be written as A/ω 2 cos(ωt + π). The phase of the excitation displacement when the position of the mass reaches a maximum (at t 1 ) thus becomes:

Phase = ωt 1 + π mod[2π] = -nψ + π mod[2π] (31)
If this phase is zero, the displacement of the excitation reaches a maximum when the position of the mass reaches a maximum. If this phase is π, the displacement of the excitation reaches a minimum when the position of the mass reaches a maximum. As the calculated phase is the phase of the excitation displacement, the result is generally positive: the frame displacement is ahead of the mass displacement.

In the model, the Fourier series of the different behaviors (x and v) were truncated to order N = 6 (to justify this choice, the reader can refer to the examples shown in Figure 5). For high orbits, subharmonics of order greater than five are not represented because they are not interesting for energy harvesting (their energy level and stability robustness are low compared to those of the behaviors plotted in Figure 9). Low orbits are not interesting for energy harvesting either. However, some of them are plotted in Figure 9 in order to know all the frequency ranges on which high orbits coexist with a low orbit. In these frequency ranges, high orbits may therefore fall on a low orbit, thus reducing the amount of energy converted by the bistable harvester.

Figure 9 shows a good agreement between the theoretical model and the experimental data. The frequency ranges of the different behaviors, their amplitude and their phases are well predicted by the theory when the stability robustness criterion is taken into account. The experimental data therefore confirm the relevance of this mathematical model which will then be exploited in the next section in order to give some guidelines for optimization of future bistable harvester prototypes.

IV. Model exploitation for optimized bistable harvesters

The previous sections have introduced a theoretical model for piezoelectric bistable harvesters able to predict their frequency response including harmonic 1 and subharmonic behaviors. This model will now be exploited to investigate the influence of different parameters (mass, stiffness, buckling level) on this frequency response to offer some leads to optimize future prototypes. A comparison with linear harvesters is also proposed all along this analysis.

The set of parameters chosen to fully describe the bistable harvester is (M, k, x 0 , L, µ, R, α, C 0 ). With this set, the general equations governing the behavior of the piezoelectric bistable harvester with low buckling level (Equation ( 5)) become:

       ẍ + x 2 0 L 2 2k M x 2 x 2 0 -1 x + µ M ẋ + 2α ML xv = -A cos(ωt) C 0 v = 2 α L x ẋ - v R (
32) In the following, not all these parameters will be considered as variable, but only the relevant ones, as explained in the followings. First, x 0 is bijectively coupled to the length L: a value of x 0 is only obtained for a single value of L and vice-versa. Therefore, the length L will not be considered here as its influence can be deduced from the influence of x 0 . The coefficient µ mostly reflects the friction of the air on the mass. This coefficient does not need to be studied: it has to be reduced as much as possible to reduce the mechanical losses in the harvester. It will therefore not be considered in the following study. Instead of α and C 0 , the piezoelectric converter will be characterized by a more global parameter, namely the electromechanical coupling factor k 2 em , representing its ability to convert the mechanical energy supplied by the source into electrical energy. It is defined in open circuit (I = 0) for a quasi-static regime (the kinetic energy of the mass is low in front of its elastic energy) by the following expression:

k 2 em =
Electrostatic energy stored in the piezo Mechanical energy supplied by the source (33)

The electromechanical coupling factor is therefore linked to two different aspects: (i) the ratio between the elastic energy stored in the piezoelectric material used in the harvester and the elastic energy stored in the other materials; (ii) the quality of this piezoelectric material (more or less able to convert its elastic energy into electrical energy). The electromechanical coupling factor is considered near-optimal for our prototype: the latter cannot be increased significantly. Indeed, (i) the piezoelectric material used is PZT (in 33 direction) which has the best ratio quality/price of nowadays market; (ii) the materials (other than piezoelectric) storing some unconverted elastic energy during the deformations are the flexible beams and the elliptical shell of the piezoelectric stacks.

The first one stores a negligible quantity of elastic energy compared to the quantity stored in the piezoelectric element. The second one stores a bigger quantity of elastic energy but this quantity is cannot be decreased without affecting the mechanical functions that this shell must satisfy (hold the piezoelectric stacks, transfer the stresses, etc.). Hence, the electromechanical coupling factor of the prototype is supposed to be optimal for the studied architecture. For bistable harvesters, the analytical expression of this electromechanical coupling factor as a function of α and C 0 is given by Equation (26).

The remaining parameters are (M, k, x 0 , R) which can all be easily changed on the prototype: M with the size of the mass, k with the ratio between the small axis and the long axis of the piezoelectric shell, x 0 with the buckling level of the harvester and R with the load plugged to the piezoelectric converter. Table 3 summarizes the set of parameters chosen to fully describe the bistable harvester, highlights those considered as changeable in this study (stars) and details the values which correspond to the prototype introduced in the experimental section. changeable parameters on the frequency response of the bistable harvester while the other parameters are kept constant. Only two behaviors are studied on these figures. These are the two most relevant behaviors for vibration energy harvesting: harmonic 1 high orbit (on which the working point maximizing the mean harvested power is located) and subharmonic 3 high orbit (offering a useful frequency range complementary to that of harmonic 1). Each of them is characterized by three quantities: the amplitude of oscillations of the mass, the stability robustness and the mean harvested power.

Despite its importance to define the behaviors of the bistable harvester, the phase of the excitation displacement when the position of the mass reaches a maximum has not been represented on these figures. The readability of the graph representing the latter is indeed not good (the curves are mixed and confused). For odd subharmonic behaviors, this phase always begins at 0 and grows monotonically until reaching π/2 for the theoretical cutting frequency (see Figure 9). The following solution was chosen in order to make this criterion appear without reducing the readability of the figures: the working points for which this phase is π/2 are represented by solid circles. These solid circles therefore correspond to the theoretical upper limit (cutting frequency) of the different behaviors.

Lastly, the frequency response of the equivalent linear harvester is also plotted in the following figures so that the performance of these two kinds of harvesters can be compared. This equivalent linear harvester has the same mass M, the same damping coefficient µ, the same excitation A cos(ωt), the same natural angular frequency ω 0 and the same electromechanical coupling factor k 2 em as the bistable harvester. Thus, the equivalent linear harvester is governed by the following behavior law:

     ẍ + ω 2 0 x + µ M ẋ + α lin M v = -A cos(ωt) C 0 v = α lin ẋ - v R (34)
with (as detailed in Equation ( 6))

ω 0 = x 0 L 4k M (35)
Moreover, the electromechanical coupling factor of a linear harvester is given by the following expression [START_REF] Badel | Nonlinear Conditioning Circuits for Piezoelectric Energy Harvesters[END_REF]:

k 2 em = α 2 lin Mω 2 0 C 0 =⇒ α lin = k 2 em Mω 2 0 C 0 (36)

IV.1. Influence of M

Figure 11 shows the influence of the mass M on the frequency response of the bistable harvester. The value of the other parameters correspond to the prototype introduced in the experimental section. They are detailed in Table 3. For each mass, the load resistance is set to ensure the impedance matching in order to maximize the maximal mean harvested power: R match = 1/(2C 0 ω 0 ) for the bistable harvester and R match = 1/(C 0 ω 0 ) for the linear harvester (for the bistable harvester the voltage of the piezoelectric converter oscillates at a frequency twice higher than the excitation as illustrated in Figure 3 and5).

When the mass M increases: (1) the amplitude of the mass oscillations for a given excitation frequency increases for both behaviors; (2) their theoretical cutting frequency (defined by a phase of π/2) also increases resulting in an increase in the theoretical bandwidth of the harvester; (3) their real cutting frequency (defined by the stability robustness limit) decreases slightly for harmonic behavior 1 and increases slightly for subharmonic behavior 3 resulting in a slight overall increase in the global bandwidth of the harvester in real conditions; (4) this is explained by the stability robustness which increases for both behaviors with the mass; (5) the mean harvested power for a given excitation frequency increases; (6) the maximum harvested power (theoretical and real) also increases.

It is then preferable to increase the mass M as much as possible (within the limits of the volume constraints imposed) in order to increase the mean harvested power of future bistable harvesters without significantly affecting its global frequency bandwidth.

IV.2. Influence of k

Figure 12 shows the influence of the stiffness k on the frequency response of the bistable harvester. The value of the other parameters correspond to the prototype introduced in the experimental section. They are detailed in Table 3. For each stiffness, the load resistance is set to ensure the impedance matching in order to maximize the maximal mean harvested power: R match = 1/(2C 0 ω 0 ) for the bistable harvester and R match = 1/(C 0 ω 0 ) for the linear harvester.

When the stiffness k increases: (1) the amplitude of the mass oscillations for a given excitation frequency decreases for both behaviors; (2) their theoretical cutting frequency (defined by a phase of π/2) increases; (3) their real cutting frequency (defined by the stability robustness limit) increases; (4) the stability robustness reaches a lower maximum value for both behaviors but reaches the robustness limit for a higher frequency; [START_REF] Wang | Dynamic and energetic characteristics of a bistable piezoelectric vibration energy harvester with an elastic magnifier[END_REF] The mean harvested power for a given excitation frequency decreases; (6) the maximum harvested power (theoretical and real) remains constant for harmonic 1 high orbit and slightly increases for subharmonic 3 high orbit.

Globally speaking, increasing the stiffness k leads to a shift of the power curves towards higher frequencies.

The stiffness can therefore be used (adjusting the ratio between the small axis and the long axis of the piezoelectric shell) to optimize the frequency range offered by the bistable harvester in such a way that it includes all the potential frequencies of the targeted vibration source. This adjustment must remain tight so that the robustness of the behaviors is not reduced too much.

IV.3. Influence of x 0

Figure 13 shows the influence of the stable position value x 0 on the frequency response of the bistable harvester. The value of the other parameters correspond to the prototype introduced in the experimental section. They are detailed in Table 3. For each value, the load resistance is set to ensure the impedance matching in order to maximize the maximal mean harvested power: R match = 1/(2C 0 ω 0 ) for the bistable harvester and R match = 1/(C 0 ω 0 ) for the linear harvester.

When the value of the stable positions x 0 increases: (1) the amplitude of the mass oscillations for a given excitation frequency increases for both behaviors; (2) the theoretical cutting frequency (defined by a phase of π/2) of harmonic 1 high orbit decreases while it remains relatively constant for subharmonic 3 high orbit; (3) their real cutting frequency (defined by the stability robustness limit) behave in the same way. Globally, the real frequency range of both behaviors decreases when x 0 increases. The frequency range of harmonic 1 high orbit is however more affected than that of subharmonic 3 high orbit. (4) The stability robustness decreases for both behaviors; (5) the mean harvested power for a given excitation frequency increases; (6) the maximum harvested power (theoretical and real) remains constant for har-Figure 11: Influence of the mass M on the frequency response of the bistable harvester for an excitation of 2.5 m.s -2 . Masses shown: 10 g (orange curves), 17.3 g (prototype value -blue curves) and 30 g (red curves). These values receptively correspond to the following natural angular frequencies: 159 rad.s -1 (25 Hz), 121 rad.s -1 (19 Hz) and 92 rad.s -1 (15 Hz). For each mass, the load resistance is set to ensure the impedance matching R match = 1/(2C 0 ω 0 ): 2.98 kΩ, 3.91 kΩ and 5.16 kΩ. For each behavior, the theoretical cutting frequency is reached when the phase of the phase of the excitation displacement is π/2 when the position of the mass reaches a maximum. The response of the equivalent linear harvester is shown in dashed curves. monic 1 high orbit and increases for subharmonic 3 high orbit.

From these observations, it can be interesting to increase the value of the stable positions x 0 (changing the buckling level of the harvester) until the overlap area between harmonic 1 and subharmonic 3 behaviors disappears. The mean harvested power would therefore be globally increased without losing too much bandwidth. For our prototype, the best configuration with this approach corresponds to the blue curves in Figure 13.

IV.4. Influence of R

The optimization of the parameter R is more specific than the other parameters studied above as it does not concern the mechanical harvester itself but the electrical circuit connected to its piezoelectric converter. Figure 14 shows the influence of R on harmonic 1 high orbit and subharmonic 3 high orbit. The value of the other parameters correspond to the prototype introduced in the experimental section. They are detailed in Table 3. The maximum mean harvested power of both behaviors considering stable orbits (solid lines + circles) is optimal when the load resistance is set to match the impedance of the piezoelectric converter at the natural angular frequency ω 0 of the harvester (red curves). This resistance is equal to R match = 1/(2C 0 ω 0 ) for bistable harvesters as the voltage of the piezoelectric converter oscillates at a frequency twice higher than the excitation (phenomenon illustrated in Figure 3 and5). When the load resistance is increased or decreased, the maximum mean harvested power of both behaviors decreases while their respective frequency ranges increases. Figure 15 shows a comparison between the frequency responses of the bistable harvester given by the orbits defined as stable and robust obtained with R match (red curves) and obtained with the resistance tuned (grey curves) to ensure the highest mean harvested power for each frequency. The resistance tuning is profitable on both harmonic 1 high orbit (maximal mean harvested power increased by 20% and frequency range increased by 14%) and on subharmonic 3 high orbit (frequency range increased by 49%). The global frequency range offered by the bistable harvester is increased by 32% with the resistance tuning.

In conclusion, the optimal resistance to maximize the mean harvested power over the entire frequency range of the bistable harvester is obtained when it is set equal to the impedance of the piezoelectric converter at the natural angular frequency ω 0 of the bistable harvester: R match = 1/(2C 0 ω 0 ) (red curves in Figure 14). Then, a resistance tuning when the excitation get beyond the real cutting frequencies (defined by the stability robustness limit) of the two behaviors will allow to increases the global frequency range offered by the harvester by roughly 30% and increase its maximal mean harvested power by roughly 20%.

IV.5. Bistable vs linear

The first observation that can be made from the comparison between linear and bistable harvesters in Figures 11, 12 and 13 concerns the mean harvested power: the maximum reached with the linear harvester (at its resonance frequency) is systematically equal to the maximum reached with the bistable generator (at the theoretical cutting frequency of the harmonic 1 high orbit Figure 13: Influence of the stable positions value x 0 on the frequency response of the bistable harvester for an excitation of 2.5 m.s -2 . Values shown: 0.50 mm (prototype value -orange curves), 0.70 mm (blue curves) and 0.90 mm (red curves). These values receptively correspond to the following natural angular frequencies: 121 rad.s -1 (19 Hz), 170 rad.s -1 (27 Hz), 218 rad.s -1 (35 Hz). For each mass, the load resistance is set to ensure the impedance matching R match = 1/(2C 0 ω 0 ): 3.91 kΩ, 2.80 kΩ and 2.18 kΩ. For each behavior, the theoretical cutting frequency is reached when the phase of the phase of the excitation displacement is π/2 when the position of the mass reaches a maximum. The response of the equivalent linear harvester is shown in dashed curves.

Figure 14:

Influence of the load resistance R on the behaviors of the bistable harvester for an excitation of 2.5 m.s -2 . Red curves correspond to the resistance ensuring the impedance matching with the piezoelectric converter for an excitation set at the natural angular frequency R match = 1/(2C 0 ω 0 ) (3913Ω). Blue curves correspond to resistances lower than R match and grey curves higher than R match . For each behavior, the theoretical cutting frequency is reached when the phase of the excitation displacement is π/2 when the position of the mass reaches a maximum.

Figure 15:

Frequency responses of the bistable harvester (excitation of 2.5 m.s -2 ) obtained for (i) a resistance ensuring the impedance matching R match = 1/(2C 0 ω 0 ) (red curves) and (ii) for a resistance tuned to get the highest mean harvested power for each frequency (grey curves). The latter is plotted only when it leads to better performances than the impedance matching. The frequency responses of the linear harvester is plotted as reference (tuned in dashed grey and for the impedance matching R = 1/(C 0 ω 0 ) in dashed red).

defined by a phase of π/2). Hence, in theory, those two harvesters have the same maximal performances.

However, the theoretical cutting frequency of the bistable harvester is not robust enough to be maintained over time in real conditions. The reachable maximum mean harvested power of the bistable harvester (robust orbits) is always lower than the maximum offered by the linear harvester. In terms of the maximum mean harvested power, the bistable harvester is therefore less efficient than its equivalent linear harvester [START_REF] Liu | A new figure of merit for wideband vibration energy harvesters[END_REF].

On the other hand, the two harvesters do not have the same bandwidth. The frequency range on which the bistable harvester recovers more than 100 µW 1 on average is systematically bigger than linear harvester's one.

The stability robustness criterion therefore balances the interest of the bistable harvester for vibration energy harvesting. Its maximum mean harvested power is greatly reduced with this criterion and becomes much lower than that of the equivalent linear harvester. However, this weakness is counterbalanced by its large bandwidth compared to that of the linear harvester. The latter remains relevant even after the inclusion of the stability robustness criterion, partly due to the exploitation of subharmonic 3 behavior. The use of one or the other of these harvesters will therefore be linked to fluctuations in the frequency of the ambient vibration source. For a source with a fixed and stable frequency, the linear harvester will be more interesting. For a source with a fluctuating frequency, the bistable harvester will be preferred. This conclusion is different from that obtained by considering the classical theoretical study of bistable harvester (without stability robustness criteria) which would have favoured the bistable harvester in all cases (source with stable or fluctuating frequency). Figures 11,[START_REF] Panyam | On approximating the effective bandwidth of bi-stable energy harvesters[END_REF][START_REF] Liu | Novel piezoelectric bistable oscillator architecture for wideband vibration energy harvesting[END_REF] and 15 also show how this comparison evolves with (M, k, x 0 , R). This evolution is summarized in Table 4 which introduces two different ratio for comparing the two harvesters: the ratio between their maximal mean harvested power and the ratio of their frequency range on which they harvest more than 100µW 1 . For the bistable harvester, these data are measured on robust and stable orbits. When M increases, the respective advantages of the two harvesters increase as well: the maximal mean harvested power of linear harvester increases relatively to bistable's one; the frequency range of bistable harvester increases relatively to linear's one. Decreasing x 0 is profitable for bistable harvester as its frequency range increase relatively to linear's ones while its maximal mean harvested power remains relatively constant. Increasing k or tuning the resistance R is profitable for bistable harvester as its maximal mean harvested power and its frequency range increase relatively 1 With this average, it is possible to power an isolated wireless sensor consuming 450 mJ per transmitted packet -22 mA at 3.3 V for 6.2 ms [START_REF] Gutiérrez | A Heterogeneous Wireless Identification Network for the Localization of Animals Based on Stochastic Movements[END_REF] -with one packet sent every 1.25 h. to linear's ones. * ratio of maximal mean harvested powers (robust bistable/linear).

• ratio of frequency ranges > 100µW 1 (robust bistable/linear).

V. Conclusion

The first part of this study introduced an analytical model for piezoelectric bistable harvester able to predict all its possible steady-state behaviors (including subharmonic behaviors) as a function of the excitation frequency. This model characterizes the different behaviors with two criteria: their stability to small disturbances and their stability robustness. The second criterion is calculated for stable high orbits and defines their ablity to handle bigger disturbances (that can be found in real conditions) without falling on a low orbit. This model has then be validated thanks to experimental data.

The second part of this study has exploited this model to investigate the influence of different parameters of the bistable harvester (M, k, x 0 ) which can all be easily changed on future prototypes: M with the size of the mass, k with the ratio between the small axis and the long axis of the piezoelectric shell, x 0 with the buckling level of the harvester. The influence of the load resistance R plugged to the piezoelectric converter has also been analyzed. Table 5 presents a summary of the results that have been obtained which can be used as a guide to design the future bistable harvesters with optimized performances for different potential applications. Considering the assumptions made in this article to get these results, these guidelines are valid for all the bistable harvesters which respect the two following constraints: (i) the behavioral equation must be a Duffing-type equation as detailed in Equation ( 5), (ii) the global electromechanical coupling must be relatively low. These two constraints are nevertheless satisfied by the majority of the bistable harveter prototypes presented so far in the literature [START_REF] Harne | A review of the recent research on vibration energy harvesting via bistable systems[END_REF].

Finally, this study focused on the comparison between bistable and linear harvesters. The stability robustness

Parameter Action

Inertial mass M Increase as much as possible within the limits of volume constraints imposed.

Stiffness k

Adjust it to fit the frequency range on the potential frequencies of the targeted vibration source (in the design introduced in the article, adjust the ratio between the small axis and the long axis of the piezoelectric shell). This adjustment must remain tight so that the robustness of the behaviors is not reduced too much.

Stable position x 0 Increase until the overlap area between harmonic 1 and subharmonic 3 behaviors disappears to increased the mean harvested power without losing too much bandwidth.

Load resistance R • Ensure the impedance matching with the piezoelectric converter (1/(2C 0 ω 0 )) to maximise the harvested power. • Tune it if the vibration source get beyond the real cutting frequencies (defined by robust orbits) for both harmonic 1 and subharmonic 3 behaviors to increase the global frequency range by roughly 30% and the maximal mean harvested power by roughly 20%.

criterion has balanced the interest of the bistable harvester for vibration energy harvesting. The maximum mean harvested power reachable by the bistable harvester is not as high as the common theory predicts and becomes much lower than the one reachable with the equivalent linear harvester. On the other hand, the bandwidth offered by the bistable harvester is much larger even after the inclusion of the stability robustness criterion (especially thanks to the use of subharmonic 3 behavior). One or the other of these harvesters will therefore be preferred depending on the vibration source. For a source presenting a high frequency fluctuation, the bistable harvester will show better results. For a source presenting a fixed and stable frequency, the linear harvester will be more relevant.

Figure 1 :

 1 Figure 1: Contribution of subharmonic behaviors and stability robustness criterion on the frequency response of electromagnetic bistable harvester [9].

Figure 2 :

 2 Figure 2: Architecture of the bistable harvester with piezoelectric coupling considered in this study. The different beams are supposed to be infinitely rigid.They are doubled on each side of the mass to prevent it from rotating (harvester with only one degree of freedom). The bistable harvester is illustrated in two different equilibrium positions: one of its two stable positions x = ±x 0 and its unstable position x = 0.

Figure 3 :

 3 Figure 3: Examples of time signal of steady-state behaviors sought (numerical study).

Figure 4 :

 4 Figure 4: Example of displacement-velocity phase portrait of the steady-state behaviors sought (numerical study) with their respective Poincaré map (black dots).

Figure 5 :

 5 Figure 5: Example of Fourier series decomposition of the steady-state behaviors sought (numerical study). On the X-axis, "f" refers to the fundamental component of the mass displacement and "h" refers to its different harmonics. The monochromatic excitation is represented by a vertical black line (displacement of the frame).

Figure 6 :

 6 Figure 6: Piezoelectric bistable harvester made with four beams buckled thanks to the nut.

Figure 7 :

 7 Figure 7: Admittance adjustment for the parameters (ω 0 , Q, C 0 , α) identification.

Figure 8 :

 8 Figure 8: Setup used to determine the experimental frequency response of the prototype. The latter is mounted on a shaker ensuring a sinusoidal acceleration of constant amplitude 2.5 m.s -2 .

Figure 9 :

 9 Figure 9: Theoretical and experimental spectra of the bistable harvester for a sinusoidal excitation of amplitude 2.5 m.s -2 : (a) amplitude of the different behaviors (b) phase of the excitation displacement when the position of the mass reaches a maximum (high orbits) (c) stability robustness (high orbits) (d) mean harvested power.

Figure 10 :

 10 Figure 10: Mass position signal (polychromatic) compared to its fundamental frequency for harmonic 1 high and low orbit (obtained numerically for a 30 Hz excitation) and for subharmonic 3 high orbit (obtained numerically for a 60 Hz excitation).

Figures 15 ,

 15 Figures 15,[START_REF] Stanton | Harmonic balance analysis of the bistable piezoelectric inertial generator[END_REF] 12 and 13 show the influence of these

Figure 12 :

 12 Figure 12:Influence of the stiffness k on the frequency response of the bistable harvester for an excitation of 2.5 m.s -2 . Stiffness shown: 79.5 kN.m -1 (orange curves), 159 kN.m -1 (prototype value -blue curves) and 318 kN.m -1 (red curves). These values receptively correspond to the following natural angular frequencies: 86 rad.s -1 (14 Hz), 121 rad.s -1 (19 Hz), 171 rad.s -1 (27 Hz). For each mass, the load resistance is set to ensure the impedance matching R match = 1/(2C 0 ω 0 ): 5.54 kΩ, 3.91 kΩ and 2.77 kΩ. For each behavior, the theoretical cutting frequency is reached when the phase of the phase of the excitation displacement is π/2 when the position of the mass reaches a maximum. The response of the equivalent linear harvester is shown in dashed curves.

Table 1 :

 1 Prototype main dimensions.

	Element	Value	Unit
	Inertial mass	17.3	g
	APX4 steel block thickness	8.0	mm
	Steel mass	16 × 16 × 8.0	mm 3
	Beams thickness	0.070	mm
	Horizontal beams length	25	mm
	Vertical beams length	17	mm
	Piezoelectric stack	20 × 5.0 × 9.0	mm
	Piezoelectric converter stiffness k	0.30	N.µm -1

Table 2 :

 2 Parameters of the bistable harvester prototype.

	Parameter	Symbol	Value	Unit
	Stable positions	±x 0	±0.50	mm
	Mass-frame distance	L	25	mm
	Inertial mass	M	17.3	g
	Natural angular frequency	ω 0	121	rad.s -1
	Mechanical quality factor	Q	87	-
	Excitation magnitude	A	2.5	m.s -2
	Excitation frequency	ω	15 -85	Hz
	Piezo converter force factor	α	0.068	N.V -1
	Piezo converter blocked capacity	C 0	1.05	µF
	Load resistance	R	7.83	kΩ
	Electromechanical coupling factor	k 2 em	2.74	%

Table 3 :

 3 Parameters describing the bistable harvester.

	Parameter	Symbol Proto	Unit
	Mass-frame distance	L	25	mm
	Damping	µ	0.023	N/m/s
	Blocked capacity Electromechanical coupling factor	C 0 k 2 em	1.05 2.74	µF %
	* Inertial mass	M	17.3	g
	* Stiffness	k	159	kN/m
	* Stable positions	±x 0	±0.50	mm
	* Load resistance	R	7s.83	kΩ

* parameters considered as easily changeable.

Table 4 :

 4 Bistable vs linear.

	Parameter	Value	Power ratio *	Freq. ratio •
	Mass	30	0.13	10.8
	M (g)	17.3	0.18	8.0
		10	0.24	5.8
	Stiffness	318	0.19	8.5
	k (kN/m)	159	0.18	8.0
		79.5	0.17	7.8
	Stable position	0.90	0.18	5.7
	x 0 (mm)	0.70	0.20	8.3
		0.50	0.18	8.0
	Resistance	R match	0.18	8.0
	R (Ω)	tuned	0.21	12.8

Table 5 :

 5 Guidelines for optimized bistable harvester.

Acknowledgment

The authors acknowledge the support of Région Auvergne-Rhône-Alpes through the ARC 4 Energies program.

VI. Appendix

The equations governing the behavior of bitable harvester detailed in Equation ( 5) is composed by a second order differential equation of the variable x and a first order differential equation of the variable v. Hence, the initial conditions needed to solve the problem are (x(0), ẋ(0), v(0)). Table 6 details the initial conditions that lead to the behaviors presented in Figures 3, 4 and5 thanks to a numerical study.

Behavior

Freq.

(Hz)


Pierre Graux
email: pierre.graux@inria.fr

Jean-François Lalande

Valérie Viet

Triem Tong

Jean-Francois Lalande
email: jean-francois.lalande@inria.fr

Obfuscated Android Application Development

Keywords: Security and privacy → Malware and its mitigation; Software reverse engineering obfuscation, mobile, application

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Obfuscation techniques help to protect an application from malicious users or competitors that want to gain knowledge about it. Understanding the internals of an application could be a legitimate need. For example, a security officer could be in charge of precisely understanding what is doing an application, before authorizing its installation in a company. On the other hand, malware developers want to protect the malicious code or some artifacts of the code (remote server URLs, ciphering keys,. . . ). Whatever the purpose, whether legitimate or not, the source code or the compiled code is manipulated to slow down or to prevent at least its static analysis.

Due to the nature of Android applications, specific obfuscation techniques using the capabilities of the development framework have been created. These techniques can work on both the source code or the compiled bytecode, in addition to the possibility to execute native libraries. Indeed, an Android application is delivered to the end user using a specific format of the Java bytecode, interpreted by a virtual machine named Dalvik. Besides, the operating system can compile parts of the bytecode before or during the execution. This flexibility makes possible to develop obfuscation techniques, working on the bytecode directly, without access to the original source code. Using combinations of these techniques makes obfuscated applications difficult to reverse.

In this paper, we first briefly review the seminal works about generating obfuscated Java programs and later, obfuscated Android applications. Then, we propose in Section 3 to pedagogically describe development techniques for helping the reader to build obfuscated applications 1 . These techniques are illustrated by snippets of code for the sake of clarity, even if the presented technique can be used directly on the bytecode. When developing new deobfuscation techniques, researchers can use these examples as unitary tests for improving their approach. Additionally, we discuss the difficulty to decide if an application has been obfuscated or not using the presented development technique. Indeed, finding artifacts of the use of a technique is often an easy task but deciding if the developer intent was to obfuscate parts of its code is a difficult problem. Finally, in Section 4 we present an overview of the use of these techniques on several application and malware datasets.

RELATED WORKS

Obfuscation techniques for Java programs have been firstly studied by Collberg et al. [START_REF] Collberg | A Taxonomy of Obfuscating Transformations[END_REF]. They distinguish the obfuscation of the control flow, the layout of the program and the manipulated data. All techniques are discussed in term of cost for the obfuscator (how it is easy to implement a transformation), of obfuscation robustness (the added time to reverse), of performance at execution time (the overhead). This seminal work helps to understand the transformations presented in this paper, such as field renaming or string encryption. Nevertheless, some transformations like reordering statements in basic blocks would be of no interest if the program is analyzed by automated tools or executed symbolically. Additionally, the developer could want to hide some assets, more than slowing down the reverse process. This explain why, when applied to Android system, developers would prefer to rely on transformations that hide the data more than complexify the flow. As Android enables to escape the Dalvik virtual machine into native code, new possibilities of obfuscation have been designed. Rastogi et al. [START_REF] Rastogi | DroidChameleon: Evaluating Android Anti-malware against Transformation Attacks[END_REF][START_REF] Rastogi | Catch Me If You Can: Evaluating Android Anti-Malware Against Transformation Attacks[END_REF] have worked on the design of new code obfuscation techniques and how to detect them. They have split the obfuscation techniques in two categories: the ones that can be defeated by static analysis tools or the ones that cannot. Only the use of reflection and bytecode encryption is presented as circumventing a static analysis as the static analyzer cannot even parse the code. Nevertheless, techniques used by packers, as shown in this paper, defeat easily the static analysis of the bytecode.

Recently, Dong et al. [START_REF] Dong | Understanding Android Obfuscation Techniques: A Large-Scale Investigation in the Wild[END_REF] contributed on the problem of detection of obfuscation techniques in Android applications. Their approach relies on two techniques: first, they parse the bytecode for collecting the used strings or identifiers (class, method and field names); second, if needed, they build on the collected data vectors of features using n-grams for preparing a classification phase. Then, they detect an obfuscation by classifying the vectors with an SVM classifier. This requires to train the classifier on labeled datasets which is a constraint we want to avoid. The obtained results on applications extracted from markets and malware repositories have been validated on manually analyzed datasets of goodware and malware. Nevertheless, the used datasets only focus on years 2016 and 2017 and do not use datasets from the literature. Finally, as they rely on a SVM classifier, it is difficult to express the reasons that make a sample classified as obfuscated or not. The detection rule, contained in the SVM classifier, is not expressed using real life elements (such as "this specific API is used"), but is a condition on vector values.

In this paper, we propose to enhance the results of Dong et al. [START_REF] Dong | Understanding Android Obfuscation Techniques: A Large-Scale Investigation in the Wild[END_REF] by proposing the following new contributions. First, we describe carefully how to create an application using development techniques that can be used for obfuscation. These techniques are listed in the left column of Table 1 and described from the most simple ones to the most advanced ones. Second, we provide static analysis methods for detecting the use of these development techniques. We either detect the use or the usage ratio if this technique can be applied on multiple parts of the code (listed in the middle column of Table 1). When possible, the analysis method is extended as a detection decision, as shown in the right column of Table 1.

BUILDING OBFUSCATED APPLICATIONS

To have a simple explanation, this section relies on a running example of code that could be found in malware sample. This code is Listing 2: String encrypted version typical of the code used by malware to communicate with its command and control server (C&C). This code is presented in Listing 1. This code has three sensitive assets: the value of two constants representing the url and the port of the C&C server and the call to the method send() of the DatagramSocket class used to communicate with the C&C server.

String encryption

3.1.1 Description. The string encryption is a technique that consists in replacing string constants by an encrypted version [START_REF] Rastogi | DroidChameleon: Evaluating Android Anti-malware against Transformation Attacks[END_REF]. Before using the constant, the code has to decrypt it. When reverse engineering an APK, all the constant strings are available2 . Thus, by encrypting them, the developer forces the analyst to understand the decryption code before accessing the original strings. Additionally, the encryption key can be hidden in the code, or downloaded dynamically, which increases reverse engineering efforts. As we can see in Listing 2, the C&C url is no more directly available because of the use of the AES256 ciphering algorithm.

Detection.

To detect the usage of string encryption, a naive approach is to compute the entropy of strings. Indeed, encryption algorithms such as AES, DES, RC4, output bytes with a high entropy close to the random distribution entropy. Thus, we propose to extract all strings of an application (excluding field, method and class names) and compute their entropy. The proposed detection method consists in counting the number of strings m with an entropy greater than an entropy threshold. Then, we say that an application has been obfuscated for n strings if m is greater than n. We discuss the choice of the threshold and we show the effect of choosing n = 1, 10 or 100 in Section 4.

Listing 3: Identifier renaming version

A more advanced detection technique have been used by Dong et al. [START_REF] Dong | Understanding Android Obfuscation Techniques: A Large-Scale Investigation in the Wild[END_REF]. It uses 3-grams to traverse the identifier names and an SVM classifier to recognize encrypted identifiers. This technique requires to train the classifier on a manually labeled dataset.

These detection techniques can be defeated by the use of encryption algorithms that does not modify the entropy, for example ROT13. Additionally, strings that are converted in integer arrays escape the detection.

Identifier renaming

3.2.1 Description. Similarly to the string encryption technique, the identifier renaming technique consists in replacing code identifiers by obfuscated ones. As the constant strings, the identifiers (package, class, method and field names) are available directly in the APK 3 . By replacing the identifiers by meaningless ones, the developer prevents the analyst to guess the purpose of an identifier. Moreover it can confuse him by using identifiers hard to distinguish or to remember, as shown in Listing 3. Some tools such as ProGuard 4 , propose to automate this renaming process.

Detection.

To detect the usage of identifier renaming, we have implemented a naive approach based on word lists. This heuristic relies on the fact that usually developers use composition of words to build their identifiers (either CamelCase or snake_case). We retrieve all the identifiers and we split them by underscore and uppercase in order to get the constitutive words of the identifier. We remove all the potential small words of size lesser than 4, that may be abbreviations. Then, we check the ratio r of the number of non-meaningful identifiers over the total, by searching them in a word list. The proposed detection method consists in comparing r with a chosen threshold. We discuss the choice of the threshold in Section 4.

Similarly to string encryption method, Dong et al. have used 3grams and an SVM classifier. Although they obtain good results [START_REF] Dong | Understanding Android Obfuscation Techniques: A Large-Scale Investigation in the Wild[END_REF], the learning phase uses applications obfuscated by Proguard and DashO, which implies that the detection method is efficient with these two tools. Using an unknown obfuscator renaming identifiers may be not detected by the learning machine. For these reasons, we prefer to rely on a detection method not based on a precise obfuscation tool.

Reflection

3.

3.1 Description. The reflection usage technique aims at hiding the methods and fields that the code is calling [START_REF] Rastogi | DroidChameleon: Evaluating Android Anti-malware against Transformation Attacks[END_REF]. Developers use

Listing 4: Reflection version

an API to access and modify at runtime internal object entities such as fields of an object. For obfuscation purpose, using reflection removes the direct references to the type or the method names that are used in the bytecode. For example in Listing 4, the methods getByName and send will not appear in the bytecode of the sendCommand method. Coupled with the string encryption technique, the use of reflection prevents the APK to be easily statically analyzed.

Detection.

Because reflection mechanisms are given by the Java API java.lang.reflect and thus cannot be renamed, its detection is trivial [START_REF] Dong | Understanding Android Obfuscation Techniques: A Large-Scale Investigation in the Wild[END_REF]: one only has to search for references to methods belonging to this package.

It has to be noted that reflection can be used for non-obfuscation purpose. Thus, it is non trivial to distinguish reflection usages that are made for obfuscation purpose and those that are not. Dong et al. [START_REF] Dong | Understanding Android Obfuscation Techniques: A Large-Scale Investigation in the Wild[END_REF] have developed a more advanced technique that analyzes instructions of the bytecode to recover the name of the class on which is invoked a reflection method, but this extra information is not sufficient to determine if the goal is to obfuscate.

Dynamic code loading

3.4.1 Description. Dynamic code loading can be used for loading dynamically some classes when not available at compile time or when building a distributed application loading remote components from a server. In addition, the bytecode can be ciphered and deciphered before loading it [START_REF] Rastogi | DroidChameleon: Evaluating Android Anti-malware against Transformation Attacks[END_REF]. DexFile and PathClassLoader provide helper classes for mapping a dex file into the memory and instanciate new objects. The updated snippet of code is given in Listing 5: after downloading a Jar file, the send() method is invoked using reflection.

Substituting totally the class loader at runtime is also possible [START_REF] Lalande | Teaching Android Mobile Security[END_REF]. Nevertheless, it requires to access and update the private attribute mClassLoader of the class LoadedApk of the application5 .

Detection.

Because dynamic code loading mechanisms are given by specific classes of the Android API, its detection is trivial [START_REF] Dong | Understanding Android Obfuscation Techniques: A Large-Scale Investigation in the Wild[END_REF] we only need to search for occurrences of all class loader classes. Once again, this detection does not distinguish between dynamic code loading used for obfuscation purpose and the others. Thus, it can be used to circumvent tools that only work on the bytecode. For exemple, an analyst that searches for the usage of network package methods, can miss the usage of libc socket, as shown in Listing 6.

Detection.

The detection of native methods is easy because they have to be declared as native in the bytecode and their implementation are located in a specific compiled library. However it is difficult to determine if native methods are used for obfuscation purpose or for a legitimate use, for example graphical routines using OpenGL ES.

Packer

3.6.1 Description. A packer is a tool than encrypts the bytecode of an APK, stores it in a resource file and removes the original bytecode. Finally, it adds an unpacking routine which is called before using the bytecode and is in charge of decrypting and loading the original bytecode. Thus, the bytecode is only available at runtime. Usually, the unpacking routine is developed using native code because the bytecode is constrained by the Dalvik virtual machine which protects its internal structure representing the loaded APK against modification. Numerous packing services 6 propose to pack Listing 7: Packer replacing the bytecode of a method application sent by customers. However, packer's internals are not documented. Thus, we give in the following more insight on how to implement a packer and we focus on the native code that unpack and replace parts of the bytecode.

As dynamic code loading can be easily detected by a dynamic tool hooking the relevant methods, packers use native code to change the Dalvik virtual machine internal structures representing the APK and the mapped bytecode of the memory. This principle can be exploited at different levels. Packers can either populate the bytecode for one method, for the whole class, or for the whole APK.

As an example, we present a packer that populate the bytecode of a method [START_REF] Lalande | Teaching Android Mobile Security[END_REF]. This application locates the packed bytecode from the file /proc/self/maps where the APK file is located in the memory. Then, it fully overwrites its bytecode. To locate the bytecode that needs to be re-written, the code navigates in internal structures of the Android runtime. Finally, to ensure to be run before any other bytecode, the unpacking routine can be implemented in a class Application and set in the Manifest under the android:name parameter of the application manifest block. The main steps of the native code that unpacks and replaces the bytecode is given in Listing 7.

Detection.

To detect the presence of packer we use APKiD [START_REF] Novella | APKiD: "PEiD" for Android Applications[END_REF], a tool that searches static artifacts of known packers in APKs. However, an "home-made" packer would not be detected. Detecting unknown packer remain an open research problem [START_REF] Liao | Automated Detection and Classification for Packed Android Applications[END_REF][START_REF] Ugarte-Pedrero | SoK: Deep packer inspection: A longitudinal study of the complexity of run-time packers[END_REF][START_REF] Xue | Adaptive unpacking of Android apps[END_REF].

OBFUSCATION USAGE IN DATASETS

In this section, we discuss the use of development techniques, used for obfuscation or not, for goodware and malware datasets.

Used datasets

The four datasets used in our experiments are: two well known but old datasets (AMD dataset [START_REF] Wei | Deep ground truth analysis of current android malware[END_REF] and the Drebin dataset [START_REF] Arp | Drebin: Effective and explainable detection of android malware in your pocket[END_REF]), one more recent (GM19 dataset [START_REF] Viet Triem | Isolating Malicious Code in Android Malware in the Wild[END_REF]) and one large (AndroZoo [START_REF] Allix | AndroZoo: Collecting Millions of Android Apps for the Research Community[END_REF]). Drebin is the oldest one: it contains 5560 samples from 2008 to 2012. AMD contains 24552 samples from 2010 to 2016. GM19 contains two balanced sets of 5000 goodware and 5000 malware with an homogeneous repartition of dates (2015-2018) and APK size. Finally, we built a custom large dataset by picking randomly 9,041 applications from AndroZoo in order to cover the largest time interval.

Packer, native, DCL, reflection

Table 2 compares the detection rates for the different techniques presented in Section 3. When a packer is detected, the tests for native method, dynamic code loading and reflection are not run because they are static techniques and would end up analyzing the unpacker code instead of the application code itself. The percentages associated to these three techniques are computed for the non-packed applications only. When comparing the results of the recent dataset MAL with the old AMD and Drebin datasets, we globally see an increase of native methods, dynamic code loading and reflection usage. This is confirmed by the experiment on the subset of Androzoo for which we represented the detection rate over the years 2008 to 2018 in Figure 1. In this figure, we also note an increase of the packer usage after 2014. The reader should note that the years 2008 to 2010 contain few samples analyzed (less than 164) which make the curves less precise for these dates. In particular, the year 2008 is based on 61 samples only and results seem biased.

When comparing goodware (GOOD dataset) and malware (MAL dataset), it is clear that native methods are more frequently used in malware samples. On the contrary, dynamic code loading and reflection does not discriminate a goodware from a malware. Indeed the usage of such development technique is not by itself an obfuscation attempt. For example, Google libraries, that are embedded in goodware APK file, contains reflection calls and class loading mechanisms even if these librairies do not use any obfuscation techniques (they are open-source).

String encryption and identifier renaming

For studying globally the identifier renaming, we represented in Figure 2, for each dataset, the percentage of APKs detected for every possible threshold. An application is said to use identifier renaming if the number of identifiers out of the dictionary divided by the total number of identifier is greater than the threshold. We observe that the curves decrease slowly until reaching a ratio threshold of 0.5. Thus, almost all applications have more than 50% of identifiers out of the dictionary which can be explained by our dictionary nature: it only contains natural language and no technical words of computer science domain. We also observe a slow down of the decreasing of the curve around the ratio threshold of 0.8. That means that less than 13% of applications have their identifiers out of our dictionary, and thus are probably obfuscated. Thus, we propose a detection threshold of 0.8 for deciding if an application has been obfuscated by renaming identifiers. 3 discusses the detection method based on the entropy for applications containing encrypted strings. We consider that an application is using string encryption if more than n strings have an entropy greater than a threshold t. For GOOD dataset, Figure 3 shows the percentage of detected applications for every threshold and for three different values of n (1, 10, 100). A threshold higher than 4.6 does not detect any application. The curves for other datasets are very similar and thus not shown in this paper. The choice for the entropy threshold and the value of n have a high impact on the number of selected applications. For example, n = 10 and threshold = 4.4 would select 20% of applications for the GOOD dataset. As we know that string encryption is not used a lot [START_REF] Dong | Understanding Android Obfuscation Techniques: A Large-Scale Investigation in the Wild[END_REF], such parameters would lead to get a high number of false positive. Moving the threshold to 4.6 results of a selection of 0% of applications. We conclude that choosing a threshold for a detection method is not reliable. Further investigations are needed with a ground truth labeled dataset.

Nevertheless, we manually checked the results obtained for string encryption and identifier renaming. Because such a verification is time consuming, we only measured that the top most scored are really true positives. Thus, we took the 5 top most detected malware of each dataset and reversed them manually. We observed 40% of true positives (60% of false positive) for string encryption and 100% of true positive (0% of false positive) for identifier renaming. The observed false positive are due to a legitimate usage of high entropy strings such as base64 encoded strings or fields containing the full english alphabet.

CONCLUSION

This paper has presented obfuscation techniques used by applications developers. We concentrated on showing how to develop an obfuscated application. Experiments show the increasing usage of these techniques. We also proposed simple detection techniques for obfuscation method that are simple to develop. They help to quantify the usage of each obfuscation technique. Nevertheless, choosing the detection thresholds that have a high true positive rate and low false negative rate remains a problem to investigate.

  Finally, we analyze both goodware and malware datasets of the literature and give a global view of development techniques over the years 2008-2018, a time period five times longer than Dong et al..

6 / * Get the dex_code_item_offset_ field of the ArtMethod class 7

 7 The hardcoded offset[START_REF] Rastogi | DroidChameleon: Evaluating Android Anti-malware against Transformation Attacks[END_REF] has been retrieved using objdump on the compiled runtime. This offset could change from one Android version to another * / 8 unsigned int code_item_offset = * (unsigned int * ) (( char * )art_method + 8); 10 / * Find the location of the DEX in the memory using / proc/self/maps. To found the DEX area, search for the .odex file * / 11 void * dex_file_location = getDexFileLocation(); 13 / * Get the new bytecode from a PNG file * / 14 const void * mmaped_file_location = GetXoredApk(); 16 / * Retrieve the CODE_ITEM structure using its offset from the beginning of the DEX. CODE_ITEM is a DEX structure and is documented * / 17 void * code_item = (void * )((char * )dex_file_location + code_item_offset); 19 / * Retrieve the bytecode size and the bytecode address by navigating in the CODE_ITEM structure * / 20 unsigned int APK_insns_size_in_code_units_ = * ( unsigned int * ) ((char * )code_item + 12); 21 void * APK_insns_ = (void * ) ((char * )code_item + 16); 23 / * Do the same for the PNG (xored APK). Recomputing code_item_offset is not needed because the offsets in the APK and in the xored APK are the same * / 24 code_item = (void * )((char * )mmaped_file_location + code_item_offset); 25 unsigned int PNG_insns_size_in_code_units_ = * ( unsigned int * ) ((char * )code_item + 12); 26 void * PNG_insns_ = (void * ) ((char * )code_item + 16); 28 / * Set the method instructions writable * / 29 void * base_addr = (void * )((char * )APK_insns_ -(( unsigned long)APK_insns_ % PAGE_SIZE)); 30 mprotect(base_addr, (size_t)((char * )APK_insns_ + APK_insns_size_in_code_units_ * 2-(char * ) base_addr), PROT_READ|PROT_WRITE|PROT_EXEC); 32 / * Un-xor and copy all new instructions. 0x42 is the xor key * / 33 unsigned int i; 34 for(i=0 ; i < APK_insns_size_in_code_units_; i++) { 35 * ((char * )APK_insns_ + 2 * i) = * ((char * )PNG_insns_ + 2 * i) ^0x42; 36 * ((char * )APK_insns_ + 2 * i + 1) = * ((char * ) PNG_insns_ + 2 * i + 1) ^0x42; 37 } 38 }

Figure 1 :

 1 Figure1: Obfuscation evolution on a subset of AndroZoo[START_REF] Allix | AndroZoo: Collecting Millions of Android Apps for the Research Community[END_REF]

  Figure 2: Identifier renaming

Figure 3 :

 3 Figure 3: Entropy for GOOD dataset

Figure

  Figure3discusses the detection method based on the entropy for applications containing encrypted strings. We consider that an application is using string encryption if more than n strings have an entropy greater than a threshold t. For GOOD dataset, Figure3shows the percentage of detected applications for every threshold and for three different values of n (1, 10, 100). A threshold higher than 4.6 does not detect any application. The curves for other datasets are very similar and thus not shown in this paper. The choice for the entropy threshold and the value of n have a high impact on the number of selected applications. For example, n = 10 and threshold = 4.4 would select 20% of applications for the GOOD dataset. As we know that string encryption is not used a lot[START_REF] Dong | Understanding Android Obfuscation Techniques: A Large-Scale Investigation in the Wild[END_REF], such parameters would lead to get a high number of false positive. Moving the threshold to 4.6 results of a selection of 0% of applications. We conclude that choosing a threshold for a detection method is not reliable. Further investigations are needed with a ground truth labeled dataset.Nevertheless, we manually checked the results obtained for string encryption and identifier renaming. Because such a verification is time consuming, we only measured that the top most

Table 1 :

 1 Overview of our analysis methods

	Development technique Analysis method	Detection
	String encryption	usage ratio	yes (threshold)
	Identifier renaming	usage ratio	yes (threshold)
	Reflection	use	no
	Dynamic code loading	use	no
	Native method	use	no
	Packer	artifacts	yes

  3 https://source.android.com/devices/tech/dalvik/dex-format 4 https://www.guardsquare.com/en/products/proguard 1 class CommandAndControl { 2 public String url = "cc.url"; 3 public int port = 4242; 4 public void sendCommand(byte[] buf, int length) {

	5	InetAddress addr = InetAddress.class
	6	.getDeclaredMethod("get"+"ByName", String.class)
	7	.invoke(null, this.url);
	8	DatagramPacket dp = new DatagramPacket(buf, length,
	9	addr, this.port);
	10	DatagramSocket.class.getDeclaredMethod("s"+"end",
	11	DatagramPacket.class)
	12	.invoke(new DatagramSocket(), dp);
	13 }}

Table 2 :

 2 Obfuscation detection for various datasets

		Total Packer Native	DCL	Reflection
	GOOD [11]	4999	3 0,06%	1266 25,58% 91,82% 4544	4735 95,68%
	MAL [11]	4991	542 10,86%	2378 57,62% 90,38% 3730	3893 94,33%
	AMD [12]	24552	31 0,13%	5206 21,52% 63,12% 15267	19184 79,31%
	Drebin [2]	5560	0 0%	1051 19,07% 26,29% 1449	3066 55,62%

Available at: https://gitlab.inria.fr/cidre-public/obfuscated-android-unit-examples

https://source.android.com/devices/tech/dalvik/dex-format

https://gist.github.com/marshall/839003

Alibaba Inc.: http://jaq.alibaba.com, Baidu: http://app.baidu.com, Bancle: https://www. bangcle.com, Ijiami: http://www.ijiami.cn, Qihoo360: http://dev.360.cn 1 void unpack_method(jclass cls, jstring methodName, jstring methodSignature) { 2 / * Get the method ID that will be overwritten * / 3 / * In fact, in Android runtime, this ID is a pointer to an internal class named ArtMethod * / 4 void * art_method = (void * ) env->GetMethodID(cls, methodName, methodSignature);