Baptiste Bergeot
email: baptiste.bergeot@insa-cvl.fr

Sergio Bellizzi

Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network

Keywords: Nonlinear energy sink network, passive mitigation, relaxation oscillations, mitigation limit, asymptotic analysis. 1

published or not. The documents may come

Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network

Introduction

A Nonlinear Energy Sink (NES) refers to an essential nonlinear absorber based on the concept of Targeted Energy Transfer (TET). It has been shown that NESs can be used successfully for vibration mitigation [START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF], noise reduction [START_REF] Bellet | Experimental study of targeted energy transfer from an acoustic system to a nonlinear membrane absorber[END_REF] or seismic mitigation [START_REF] Nucera | Targeted energy transfers in vibro-impact oscillators for seismic mitigation[END_REF]. Recently, particular attention has been paid to the comparison between a single NES versus several NESs. In [START_REF] Li | Dynamics of two vibro-impact nonlinear energy sinks in parallel under periodic and transient excitations[END_REF], a linear oscillator coupled with two vibro-impact NESs in parallel is studied under periodic and transient excitations, respectively. Several coupled parallel NES systems involving cubic nonlinearities are considered in [START_REF] Vaurigaud | Targeted energy transfer with parallel nonlinear energy sinks. Part I: Design theory and numerical results[END_REF] to control strong modes of a linear master structure under periodic and transient excitations, respectively. Associated experimental results are reported in [START_REF] Alireza | Targeted energy transfer with parallel nonlinear energy sinks, part II: Theory and experiments[END_REF]. Taghipour and Dardel [START_REF] Taghipour | Steady state dynamics and robustness of a harmonically excited essentially nonlinear oscillator coupled with a two-DOF nonlinear energy sink[END_REF] compare the steady state dynamical behavior of a two Degrees Of Freedom (DOFs) system -composed of an harmonically excited nonlinear oscillator coupled with a single DOF NES -with the behavior of a system consisting of a nonlinear oscillator coupled with a two-DOFs NES subjected to external harmonic excitations. Conclusions of this paper illustrate that increasing the DOFs of the NES would increase the robustness of the system to the changes in system parameters and the amplitude of the external forces.

NESs can also be used to mitigate dynamic instabilities. In [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF], the self-excitation response regimes of a Van der Pol oscillator coupled to a NES are investigated. Periodic responses, global bifurcations of different types and basins of attraction of various selfexcitation regimes are revealed using an asymptotic analysis of the system related to the decomposition of the averaged flow into slow/super-slow time scales. The case of a Van der Pol-Duffing system is considered in [START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in Van der Pol-Duffing oscillator with nonlinear energy sink[END_REF]. A series of papers [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF][START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 2: Experiments[END_REF][START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF] demonstrated that a NES coupled to a rigid wing in subsonic flow can partially or even completely suppress aeroelastic instabilities. Aeroelastic instability suppression is also considered in [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF] based on a general nonlinear multi DOFs system. More recently, other types of instabilities have been considered. In [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF][START_REF] Bergeot | Passive suppression of helicopter ground resonance using nonlinear energy sinks attached on the helicopter blades[END_REF], a theoretical/numerical analysis of the capacity of a NES to control helicopter ground resonance instability has been performed whereas the problem of passive control of friction-induced vibrations due to mode coupling instability in braking systems using two NESs is studied in [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF].

More generally, the discussion on the relationship between the dimensionality of the super-slow manifold (hereafter referred as critical manifold), the structure of the fixed points and the observed response regimes is explored in the review paper by [START_REF] Gendelman | Targeted energy transfer in systems with external and self-excitation[END_REF]. Some of the previous cited works [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF][START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in Van der Pol-Duffing oscillator with nonlinear energy sink[END_REF][START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF] provide theoretical criteria for the prediction of the response regimes resulting from the coupling of a primary dynamical system -potentially unstable -with one NES. In all these works, the system under study is (or is reduced to) of a one-DOF primary system coupled to one NES. This leads to a one dimensional critical manifold which is essential to make the theoretical prediction relatively easy. Otherwise, if the system has a critical manifold with dimension strictly higher than one, the prediction becomes very hard to obtain. This point has been already discussed in previous works by the authors [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF][START_REF] Bergeot | Passive suppression of helicopter ground resonance using nonlinear energy sinks attached on the helicopter blades[END_REF][START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF] in the context of dynamic instability mitigation and also when a NES is used as a passive control device for a nonlinear elastic string, in internal resonance conditions, excited by an external harmonic force [START_REF] Luongo | Nonlinear energy sink to control elastic strings: the internal resonance case[END_REF].

The present work studied a one-DOF nonlinear unstable primary system, which undergoes limit cycle oscillations with high amplitudes (hereafter referred to as harmful situations), coupled to several parallel NESs. The asymptotic analysis of the system shows that the critical manifold of the system can be reduced to a one dimensional parametric curve. Therefore, the analytical prediction of the response regimes can be performed. In particular, we introduce the mitigation limit of the NESs, defined as the value of a chosen bifurcation parameter which separates harmful situations (corresponding to responses in which the NESs do not act) from harmless situations (corresponding to responses in which the NESs act). The mitigation limit is therefore predicted and a literal expression, function of the NESs parameters, is even obtained in the case of identical NESs and a linearized primary system.

The next section introduces the system under study. In Sect. 3, the asymptotic analysis of the system, which associates complexification-averaging method and geometric singular perturbation theory [START_REF] Christopherk | Geometric singular perturbation theory[END_REF], is performed in the general case of a nonlinear primary system coupled to N different NESs. The particular case of a linearized primary system coupled to N identical NESs is analyzed in Sect. 4. Sect. 5 is devoted to illustrating and validating analytical predictions proposed in Sects. 3 and 4. For this purpose, a Van der Pol oscillator is considered as a primary system.

System under study

One considers a one degree of freedom unstable nonlinear system

m d 2 x dt 2 -p dx dt + kx + fNL x, dx dt = 0, ( 1
)
where m is the mass, k is the linear stiffness and -p (with p > 0) is the negative damping which causes the linear instability. The nonlinear function fNL (called nonlinearity of the primary system) allows the existence of Limit Cycle Oscillations (LCOs) on which the system can saturate. In the remainder of the paper, we assume large amplitude LCOs which correspond to harmful situations for the system.

A network of N NESs with masses m n , damping coefficients c n and pure cubic stiffnesses k N L n (for n = 1, . . . , N ), are attached to the system in an ungrounded configuration (NESs are parallel to each other, see Fig. 1). Because a NES is an essentially nonlinear oscillator, the linear stiffness components (k n ) are neglected. This assumption is in agreement with experimental data (see for example [START_REF] Bellet | Experimental study of targeted energy transfer from an acoustic system to a nonlinear membrane absorber[END_REF][START_REF] Gourdon | Nonlinear energy pumping under transient forcing with strongly nonlinear coupling: Theoretical and experimental results[END_REF][START_REF] Gaëtan Kerschen | Experimental demonstration of transient resonance capture in a system of two coupled oscillators with essential stiffness nonlinearity[END_REF]).

Taking into account the NESs displacements h n (t ), the equations of motion of the coupled system become

m d 2 x dt 2 -p dx dt + kx + fNL x, dx dt + (2a) N n=1 c n dx dt - dh n dt + k N L n (x -h n ) 3 = 0 m n d 2 h n dt 2 -c n dx dt - dh n dt -k N L n (x -h n ) 3 = 0, (2b)
for n = 1, . . . , N.

Introducing the new time t = ωt , with ω = k/m, and rescaling system of Eqs. (2) through a dimensionless small parameter > 0 assuming that the mass of the NESs are small with respect to the mass of the primary system, one obtains the following rescaled system where

ẍ -ρ ẋ + x + f N L (x, ẋ) + (3a) N n=1 µ n ẋ -ḣn + α n (x -h n ) 3 = 0 a n ḧn -µ n ẋ -ḣn -α n (x -h n ) 3 = 0, ( 3b
) for n = 1, . . . , N, m x(t) k -p f N L (x, ẋ) NES 1 NES n NES N . . . . . . (a) NES n c n k N L n h n (t) m n (b)
s ρ = p √ m k , f N L (x, ẋ) = fNL (x, ω ẋ) mω 2 , a n = m n m , µ n = c n √ m k , α n = k N L n mω 2 , ( 4
) and the dot denotes the differentiation with respect to the new time variable t.

3 Asymptotic analysis in the general case of a nonlinear primary system coupled to N different NESs

The slow-flow

First, to simplify the following calculations, it is convenient to introduce new coordinates as

u = x + N n=1 a n h n (5a) v n = x -h n , for n = 1, . . . , N, (5b)
giving reciprocally,

x = u + N n=1 a n v n 1 + N n=1 a n (6a) h n = u -1 + N k=1 k =n a k v n + N k=1 k =n a k v k 1 + N k=1 a k , (6b)
for n = 1, . . . , N.

Using Eqs. ( 5) and ( 6), system of Eqs ( 3) is transformed to the following form

ü + u + -ρ u + N n=1 a n (v n -u) + (7a) f N L (u, u, v 1 , v1 , . . . , v N , vN ) = 0 vn + u + µ n a n vn + α n a n v 3 n + (7b) -ρ u + N k=1 a k (v k -u) + µ k vk + α k v 3 k + g n,N L (u, u, v 1 , v1 , . . . , v N , vN ) = 0, for n = 1, . . . , N
after performing a first-order Taylor expansion around = 0. Because of the change of variables ( 5), the nonlinearity of the primary system affects both variables u and v n . To simplify the notations, we use again f N L even if the function has changed from Eq. (3a) to Eq. (7a).

We restrict our attention to the motion of the system in the vicinity of a 1:1 resonance associated to a frequency close to the natural frequency (ω = 1) of the primary system.

Following the Complexification-Averaging method (CA-X) [START_REF] Manevitch | Complex representation of dynamics of coupled nonlinear oscillators[END_REF][START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF], complex variables are introduced as

ψ = u + ju (8a) ζ n = vn + jv n , for n = 1, . . . , N (8b)
and

ψ = φe jt (9a) ζ n = ξ n e jt , for n = 1, . . . , N (9b)
where j 2 = -1 and φ and ξ n (for n = 1, . . . , N ) are the complex slow modulated amplitudes of the fast component e jt . Substituting Eqs. ( 8) and ( 9) into [START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in Van der Pol-Duffing oscillator with nonlinear energy sink[END_REF] and averaging over one period, system of Eqs [START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in Van der Pol-Duffing oscillator with nonlinear energy sink[END_REF]

reduces to φ = f (φ, ξ 1 , . . . , ξ N ) (10a) ξ n = g n (φ, ξ 1 , . . . , ξ N , ) , for n = 1, . . . , N, (10b)
where

f (φ, ξ 1 , . . . , ξ N ) = 1 2 ρ φ + j N n=1 a n (ξ n -φ) + f N L (φ, ξ 1 , . . . , ξ N ) , ( 11a
)
g n (φ, ξ 1 , . . . , ξ N , ) = j 2 φ - 1 2 µ n a n + j ξ n + 3jα n ξ n |ξ n | 2 8a n + 8 4 ρ -j N k=1 a k φ- N k=1 4 (µ k -ja n ) ξ k -3jα k ξ k |ξ k | 2 + g n,N L (φ, ξ 1 , . . . , ξ N ) . (11b)
Without loss of generality, the same notations have been used to express the nonlinear functions f N L and g n,N L before and after the averaging process. system of Eqs [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] describes the Complex Form of the Slow-Flow (CFSF) of [START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in Van der Pol-Duffing oscillator with nonlinear energy sink[END_REF]. It is important to note that due to the assumption of 1:1 resonance, only the terms with fast frequency equal to 1 are kept during the averaging step. Therefore, only odd nonlinearities can be taken into account, the even nonlinearities are averaged out.

The Critical Manifold

The slow-flow described by system of Eqs [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] highlights explicitly the "slow/fast" nature of the system: Eq (10a) contains only O( ) terms and Eq (10b) contains both O(1) and O( ) terms. Note that the nonlinearities of the primary system in Eq (10b) are located only in O( ) terms (see Eq. (11b)). Here the terminology used by Gendelman et al. (see e.g. [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF]) is preferred, i.e. the terms fast and slow are replaced by slow and super-slow respectively, whereas the term fast denotes the time scale determined by fast oscillations of the primary system with frequency 1.

Eq. ( 10) can be reformulated by switching from the slow time scale t to the super-slow time scale τ = t as

φ = f (φ, ξ 1 , . . . , ξ N ) (12a) ξ n = g n (φ, ξ 1 , . . . , ξ N , ) , for n = 1, . . . , N, (12b)
where = d dτ . Solutions of the super-slow/slow system (10) (or [START_REF] Gourdon | Nonlinear energy pumping under transient forcing with strongly nonlinear coupling: Theoretical and experimental results[END_REF]) can exhibit slow and super-slow epochs characterized by the speed at which the solution advances.

Stating = 0, the following subsystems are derived from [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] and [START_REF] Gourdon | Nonlinear energy pumping under transient forcing with strongly nonlinear coupling: Theoretical and experimental results[END_REF]

respectively φ = 0 (13a) ξn = g n (φ, ξ 1 , . . . , ξ N , 0) , for n = 1, . . . , N, (13b)
which is the slow subsystem, and

φ = f (φ, ξ 1 , . . . , ξ N ) (14a) 0 = g n (φ, ξ 1 , . . . , ξ N , 0) , for n = 1, . . . , N, (14b)
which is the super-slow subsystem.

The algebraic Eqs. (14b) of the super-slow subsystem defines the so-called Critical Manifold S [START_REF] Christopherk | Geometric singular perturbation theory[END_REF] as

S := (φ, ξ 1 , . . . , ξ N ) ∈ C N +1
g n (φ, ξ 1 , . . . , ξ N , 0) = 0 for n = 1, . . . , N . [START_REF] Gaëtan Kerschen | Experimental demonstration of transient resonance capture in a system of two coupled oscillators with essential stiffness nonlinearity[END_REF] Looking at Eq. (11b), one can see that Eq. (14b) can take the following equivalent form

φ = ξ n F n (|ξ n |) , for n = 1, . . . , N, (16)
where the complex functions F n are defined by [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF] showing that the form of S depends, at the time t = ωt ,only on the NESs parameters. It is convenient to characterize S in the real domain. For this purpose, polar coordinates are introduced φ = se jδ (18a)

F n (X) = 1 - 3α n 4a n X 2 -j µ n a n = R n (X) + jI n (X)
ξ n = r n e jθn , for n = 1, . . . , N, (18b)
and the modulus and the argument of Eq. ( 16) are successively computed, leading to the following equivalent form

s 2 = H n (r n ) (19a) ϑ n = -arg (F n (r n )) for n = 1, . . . , N (19b)
where ϑ n = θ nδ denotes the difference of the angular coordinate between ξ n and φ and

H n (r n ) = r 2 n R n (r n ) 2 + I n (r n ) 2 . ( 20
)
Result 3.1 Due to system of Eqs [START_REF] Li | Dynamics of two vibro-impact nonlinear energy sinks in parallel under periodic and transient excitations[END_REF], the critical manifold S is characterized as a one dimensional parametric curve evolving in R + N . It is the projection on the (r 1 , . . . , r N )-space (i.e. R + N ) of the intersection of the N Eqs. (19a), s 2 = H 1 (r 1 ), . . . , s 2 = H N (r N ), defined in the (s, r 1 , . . . , r N )-space (i.e. R + N +1 ).

Note that in the sequel we will abundantly use the following relationships holding in S

ξ n F n (|ξ n |) = ξ m F m (|ξ m |) , for n, m = 1, . . . , N (21)
or, using polar coordinates [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 2: Experiments[END_REF],

H n (r n ) = H m (r m ), θ nm = θ n -θ m = arg F m (r m ) F n (r n ) ,
for n, m = 1, . . . , N. [START_REF] Manevitch | Complex representation of dynamics of coupled nonlinear oscillators[END_REF] An illustration of Result 3.1 is presented in Fig. 2 for N = 2. One can see a typical critical manifold plotted in the (s, r 1 , r 2 )-space (i.e. R + 3 ) and its projection on the (r 1 , r 2 )-plane (i.e. R + 2 ).

Each point of S is a fixed point for the slow subsystem [START_REF] Grasman | Asymptotic Methods for Relaxation Oscillations and Applications[END_REF]. To determine the stability of a point of S, one must know if it attracts or repels the slow dynamics. For this purpose, using again the polar coordinates ( 18) combined with Eqs. ( 16), ( 17), ( 19) and ( 20), system of Eqs ( 13) can be rewritten as follows

ṡ = 0 (23a) ṙn = 1 2 (s sin ϑ n + r n I n (r n )) (23b) θn = 1 2 s r n cos ϑ n -R n (r n ) , for n = 1, . . . , N. ( 23c
)
Stability range is then determined by examining the sign of the eigenvalue real parts of the Jacobian matrix of the differential system (23b)-(23c) on critical manifold S. It can be shown that the condition of stability for a point r = [s, r 1 , . . . , r N ] is equivalent to

dH n (r n ) dr n > 0, for n = 1, . . . , N. ( 24
)
For each n, the local extrema of the real functions H n are given by the positive roots of its derivative dH n (r n )/dr n . An easy calculus shows that the local extrema occur at

r M n = 2 3 2a n -a 2 n -3µ 2 n α n (25a) r m n = 2 3 2a n + a 2 n -3µ 2 n α n , ( 25b
)
if the following relation holds

µ n < a n √ 3 , ( 26
)
and in this case r M n < r m n . Assuming that the condition (26) holds for the N NESs, the stability domain on S is only characterized by the (r 1 , r 2 , . . . , r n )-subspace D given in terms of Cartesian product as

D = N n=1 D n , ( 27
)
where

D n = 0 r M n ∪ [r m n + ∞) . ( 28
)
D n denotes the stability range associated to the n-th NES characterized on R + by the two roots, r m n and r M n (see Eqs. (25b) and (25a)).

A typical projection of the Critical Manifold onto the subspace R + 2 corresponding to the n-th NES (variables (s, r n )) is depicted in Fig. 3 where the stability range is also reported. The two points (s m , r m n ) and (s M , r M n ) where the two scalars

s m n = H n (r m n ) and s M n = H n (r M n )
characterize the bounds where S ceases to be hyperbolic 1 connecting stable or attractive (continuous line) and unstable or repulsive (dashed line) parts of S. These two points are called fold points. The two scalars r d n and r u n which will later be used are solutions of

H n (r m n ) = H n r d n ⇒ r d n = 2 √ 2 3 a n -a 2 n -3µ 2 n α n , ( 29a
)
H n r M n = H n (r u n ) ⇒ r u n = 2 √ 2 3 a n + a 2 n -3µ 2 n α n . ( 29b
)
Eq. ( 27) may be expanded as the union of 2 N Cartesian products, each of them being the product of close or left-close intervals resulting in an isolated stable or attractive part of the the critical manifold S, hereafter denoted

I a k (with k ∈ [1, 2 N ]) D = 2 N k=1 I a k , ( 30
)
where the first and the last terms are chosen as

I a 1 = N n=1 [0, r M n ] and I a 2 N = N n=1 [r m n , +∞) respectively.
Examples of typical critical manifolds for a system containing two and three parallel NESs are presented in Figs. 4(a) and 4(b) respectively. In both figures the stability domain of S is shown. In particular, Eq. ( 30) 1 S is hyperbolic if all the eigenvalues of the Jacobian matrices of the differential system [START_REF] Nucera | Targeted energy transfers in vibro-impact oscillators for seismic mitigation[END_REF] have nonzero real parts.

is illustrated in Fig. 4(a), in the (r 1 , r 2 )-plane. Indeed, D, the stable (or attractive) domain (gray areas on the figure in which S is depicted in magenta) is defined by

D = 0 r M 1 ∪ [r m 1 + ∞) × 0 r M 2 ∪ [r m 2 + ∞) = 0 r M 1 × 0 r M 2 ∪ 0 r M 1 × [r m 2 + ∞) ∪ [r m 1 + ∞) × 0 r M 2 ∪ ([r m 1 + ∞) × [r m 2 + ∞)) , ( 31
)
whereas the unstable or repulsive part (white areas on the figure in which S is depicted in blue) is defined as the associated complement set of

D in R + 2 .
Such structures of the critical manifold S provide a possibility for relaxation oscillations [START_REF] Grasman | Asymptotic Methods for Relaxation Oscillations and Applications[END_REF] of the slowflow characterized by sudden transitions (jumps) of the dynamics during each cycle (the hypothetical sudden transitions between the two stable branches are denoted by arrows in Fig. 3). Unlike the case with just one NES, the presence of several NESs allows the existence of several relaxation oscillations scenarii: relaxation oscillations in all planes (r n , s) or just in some of them. Such relaxation oscillations of the slow-flow reflect the presence of TET from the primary unstable system to the NESs and explain the existence of Strongly Modulated Responses [START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF][START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF][START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] (SMRs) of the non-averaged system [START_REF] Benoît | Chasse au canard ("duck hunting[END_REF]. In view of the foregoing, in order to enable the existence of such SMRs, the NES must respect the condition [START_REF] Alireza | Targeted energy transfer with parallel nonlinear energy sinks, part II: Theory and experiments[END_REF] to ensure the S-shape of the critical manifold S.

We proceed to a more detailed analysis of the possible steady-state regimes in the next section computing fixed points of the slow-flow and their stability.

The fixed points of the slow-flow

We are interested in the nontrivial fixed points of the slow-flow [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF]. These nontrivial fixed points characterize periodic solutions of system of Eqs [START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in Van der Pol-Duffing oscillator with nonlinear energy sink[END_REF] having a frequency exactly equal to 1, the frequency used to define the complex variables [START_REF] Gendelman | Targeted energy transfer in systems with external and self-excitation[END_REF]. The fixed point of [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] can be computed from the Real Form of the Slow-Flow (denoted RFSF) obtained by substituting the polar coordinates [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 2: Experiments[END_REF] and considering instead of the arguments δ and θ n (t), the argument differences ϑ n = θ nδ (the master component can be chosen arbitrarily, φ is here chosen for convenience). RFSF may be formally written as follows

Ż = G (Z, ) with Z = [s r 1 . . . r N ϑ 1 . . . ϑ N ]
t . (32) system of Eqs (32) gives a formulation to characterize the existence of stable periodic response which is not easy to manipulate and solve.

In the sequel, the Geometric Singular Perturbation Theory (GSPT) [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF][START_REF] Christopherk | Geometric singular perturbation theory[END_REF][START_REF] Desroches | Mixed-mode oscillations with multiple time scales[END_REF] is preferred to conduct the analysis. The dynamics of the slow-flow (CFSF [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] or RFSF (32)) for 0 < 1 is described from the analysis of the slow and super-slow subsystems (which are defined for = 0). More precisely, we use the following result of the GSPT: if 0 < 1, the dynamics of the slow-flow during the slow (resp. super-slow) epoch is given by the dynamics of the slow (resp. super-slow) subsystem (13) (resp. ( 14)). This allows to easily compute the fixed points of system of Eqs (32) which are assumed to exist only at the super-slow time scale.

As detailed in Appendix A, Eq. (32) gives access to the slow subsystem (system of Eqs (64)) and to the super-slow subsystem (system of Eqs (65)) in polor coordinates. Substiting system of Eqs [START_REF] Li | Dynamics of two vibro-impact nonlinear energy sinks in parallel under periodic and transient excitations[END_REF] with n = 1 (any n ∈ [1, N ] can be chosen) into the super-slow subsystem we obtain

dH 1 (r 1 ) dr 1 r 1 = f r1 (r 1 , . . . , r N ) , ( 33
)
where H 1 is defined by Eq. ( 20) and f r1 is built in Appendix A and explicitly given here after. Eq. ( 33) describes the super-slow dynamics projected on the critical manifold S.

From Eq. (33), it is possible to detect fixed points as

f r1 (r 1 , . . . , r N ) = 0 (34a) dH 1 (r 1 ) dr 1 = 0, (34b)
and folded singularities as

f r1 (r 1 , . . . , r N ) = 0 (35a) dH 1 (r 1 ) dr 1 = 0. ( 35b
)
Folded singularities correspond to situations for which fixed points and fold points coincide, they are hints of possible canard explosions [START_REF] Benoît | Chasse au canard ("duck hunting[END_REF]. The rest of the section is dedicated to the analysis of the fixed points which, for 0 < 1, are assumed to be the fixed points of (32).

Eq. ( 33) can be re-written separating the linear part of the primary system and the nonlinearities due to the NESs attachments as

f r1 (r 1 , . . . , r N ) = 2 ρH 1 (r 1 )- N k=1 r 2 k µ k +f N L (r 1 , . . . , r N ) , (36)
where, without loss of generality, the same notations have been used to express the nonlinear function f N L (of the primary system) before and after the use of polar coordinates [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 2: Experiments[END_REF].

Hence, finding the fixed points of system of Eqs (32) amounts to finding the roots of

ρH 1 (r 1 ) - N k=1 r 2 k µ k + f N L (r 1 , . . . , r N ) = 0 (37a) H n (r n ) -H 1 (r 1 ) = 0, for n = 2, . . . , N (37b)
which can be easily solved with a regular computer (at least up to N ≈ 5 or 6).

To check the stability of a fixed point (denoted r * = [r * 1 , . . . , r * N ]), Eq. ( 33) is written as

r 1 = f (r 1 , . . . , r N ) , ( 38
)
where

f (r 1 , . . . , r N ) = 2 ρH 1 (r 1 ) - N k=1 r 2 k µ k + f N L (r 1 , . . . , r N ) dH1(r1) dr1
.

(39) The stability is deduced from the sign of df dr1 (r * ) where, using Eq. (37b), df dr1 (r) takes the form

df dr 1 (r) = ∂f ∂r 1 (r) + N n=2 ∂f ∂r n (r) dr n dr 1 = ∂f ∂r 1 (r) + N n=2 ∂f ∂r n (r) H 1 (r 1 ) H n (r n ) . ( 40
)
Result 3.2 A fixed point r * of Eq. (32) is stable if the two following conditions are satisfied

1. df dr 1 (r * ) < 0, 2. ∀n ∈ [1, N ], r * n ∈ [0 r M n ] ∪ [r m n + ∞). A stable (respectively unstable) fixed point is de- noted r * s = r * s,1 , . . . , r * s,N (respectively r * u = r * u,1 , . . . , r * u,N
). The second condition means that to be stable the fixed point must be on a stable branch of the critical manifold S (see Sect. 3.2).

Prediction of the mitigation limit

Four main types of response regimes may be generated when a NES is attached to an unstable system: complete suppression of the instability, mitigation through Periodic Response (PR), mitigation through Strongly Modulated Response (SMR) or no mitigation. These four types of responses have been observed and analyzed in the case of a Van de Pol oscillator coupled to one NES [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF]. They have also been observed in [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF] and studied theoretically in [START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF] in the context of the mitigation of aeroelastic instabilities of a rigid wing in a subsonic flow. Furthermore, Bergeot et al. [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF] observed these responses studying mitigation of a mode-coupling instability in breaking systems using two NESs.

Numerical simulations by the authors on the considered system show that the same regimes are encountered considering N parallel NESs. In this study we classify these regimes into two categories depending on whether the NESs act or not and therefore separating harmless situations from harmful situations:

1. Harmless situation: The NESs act, resulting in the following regimes (a) Complete suppression. The trivial fixed point of the slow-flow is reached. In this case, due to the additional damping, the NESs attachments stabilize the system, i.e. the trivial fixed point, common to both the nonaveraged system and the slow-flow, becomes stable.

(b) Mitigation through Periodic Response. A nontrivial stable fixed point of the slow-flow is reached. In this case, the steady-state regime of the non-averaged system (2) is a periodic regime, i.e. LCOs with amplitudes smaller than the amplitudes of the LCOs undergone by the primary system alone.

(c) Mitigation through Strongly Modulated Response. In this case, the steady-state regime is a quasiperiodic regime which exhibits a "fast" component and a "slow" component corresponding to the envelope of the signal. The term "Strongly modulated response" has been introduced in [START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF] for the study of a harmonically forced linear system coupled to a NES. SMRs for (2) corresponds to relaxation oscillations of the slow-flow mentioned in Sect. 3.2.

To ensure the existence of SMRs, at least one of the N parallel NESs must respect the condition (26).

Harmful situation:

The NESs do not act resulting in (a) No mitigation. The NESs are not able to mitigate the instability and the non-averaged system saturates on a LCO which has an amplitude close to that of the case without NES.

Regarding the slow-flow, it reaches a stable fixed point with large amplitude.

The nature of the steady-state regime depends on two characteristics of the slow-flow: (i) the initial conditions (to know where the slow dynamics leads the trajectory on S) and (ii) the fixed points (position and stability) which orientate the dynamics on S at the super-slow time scale. In the present paper, for convenience and in accordance with real word situations, we consider a set of initial conditions [s(0), r 1 (0), . . . , r N (0)] as a small perturbation of the trivial solution.

In reference to the two situations (harmless and harmful) described above and the previous comments, the mitigation limit is introduced in the following definition.

Definition 3.1 Considering a set of initial conditions (for the slow-flow)

[s(0), r 1 (0), . . . , r N (0)] as a small perturbation of the trivial solution, the mitigation limit is defined as the value of a chosen bifurcation parameter which separates harmful situation from harmless situation.

Classically the bifurcation parameter is one of the parameters of the considered system which affects the fixed points of the slow subsystem (see system of Eqs (37)).

The mitigation limit is only based on the nature of the steady-state regimes when 1:1 resonance occurs. This approach may be questioned regarding the amplitudes of these regimes. In this context, two issues have to be given a special attention. First, one might wonder, if a stable fixed point of the slow-flow is reached, how to know if the system undergoes a Mitigation through Periodic Response (harmless situation) or No mitigation (harmful situation). In general, the primary nonlinearity generates, in the bifurcation diagram of the slow-flow, an isolated branch of stable fixed points which are the larger ones (compared to both stable and unstable fixed points). In particular situations (addressed in Sect. 5.3) some other stable fixed points may also have large amplitude, remaining smaller than those of the primary nonlinearity. Consequently, the following issue may arise: if SMRs are observed, is the resulting amplitude sufficiently small to characterize the regime as a harmless situation? This point is easier to answer. Indeed, the amplitude of the relaxation oscillation of the slow-flow is given by the coordinates of the fold points (29). Therefore, the NESs parameters can be chosen to have relaxation oscillations with small amplitudes compared with the amplitudes of the LCOs generated by the primary nonlinearity.

In the following, the method to predict the steadystate regimes, and consequently the mitigation limit, is presented. It is based on the above study (Sects. 3.2 and 3.3) which provides a theoretical framework to explain and predict the possible steady-state regimes (described above) of first the slow-flow itself and consequently of the non-averaged system [START_REF] Benoît | Chasse au canard ("duck hunting[END_REF].

Conditions to ensure that the system undergoes a harmless situation are stated in the following result.

Result 3.3 The conditions to obtain harmless situations are split in two cases:

Case 1 We assume that at least one of the N parallel

NESs respects condition [START_REF] Alireza | Targeted energy transfer with parallel nonlinear energy sinks, part II: Theory and experiments[END_REF].

Case 1a:

We assume that the slow-flow has a stable fixed point r * s = r * s,1 , . . . , r * s,N (trivial or nontrivial) on the first attractive part of the critical manifold S, i.e.

r * s ∈ I a 1 . ( 41
)
The system is in a harmless situation (trajectories reach inevitably r * s ) and undergoes Complete suppression for the trivial fixed point or Mitigation through Periodic Response for the nontrivial fixed points. If (41) is not respected, Case 1b holds instead of Cases 1a.

Case 1b: Let R u be the set of all nontrivial unstable fixed points on a stable part of the critical manifold S. We assume that R u = ∅, i.e the slow-flow has at least one nontrivial unstable fixed point r * u = r * u,1 , . . . , r * u,N on a stable part of the critical manifold S.

Let N = {n ∈ [1, N ] such that µ n < a n / √ 3}
(here N = ∅). The system is in an harmless situation if

∀n ∈ N , max r * u ∈Ru r * u,n > r u n , ( 42
)
where r u n is defined by (29b).

In this case, harmless situations correspond to Mitigation through Periodic Response or Mitigation through Strongly Modulated Response.

To discriminate PR from SMR, the existence of the stable fixed point is not sufficient. Indeed, the relaxation oscillations may prevent the trajectory to reach the latter.

Conversely, an harmful situation is obtained if the system is neither in Case 1a, i.e. condition (41) is not respected, nor in Case 1b, i.e.

• the system has no unstable fixed point. In this situation, the system reaches the stable fixed point with large amplitude and it is therefore in an harmful situation.

• or ∃n ∈ N , max

r * u ∈Ru r * u,n < r u n . ( 43
)
In this case, after a transient regime at the slow time scale the trajectory is pushed to a stable fixed point with a large amplitude (due to the primary nonlinearity and therefore close to that of the case without NES), along S (at the super-slow time scale), by the unstable fixed point r * u (which satisfies (43)) and the system is actually in an harmful situation.

Case 2 We assume that

∀n ∈ [1, N ], µ n > a n / √ 3.

In this case relaxation oscillations cannot happen. A harmful situation is obtained if the slow-flow admits only one stable nontrivial fixed point. It corresponds to the fixed point due to the primary nonlinearity. A harmless situation is obtained if at least two stable nontrivial fixed points exist.

The trajectory reaches inevitably the first one on S. Therefore, the system undergoes Complete suppression for the trivial fixed point or Mitigation through Periodic Response for the nontrivial fixed point.

A quantitative characterization of the mitigation limit of the system is now given. As an example, we consider the parameter ρ (see ( 3)) as the bifurcation parameter. Let ρ lim be the value of the parameter ρ after which the slow-flow has no more nontrivial unstable fixed point.

Result 3. [START_REF] Bergeot | Passive suppression of helicopter ground resonance using nonlinear energy sinks attached on the helicopter blades[END_REF] The mitigation limit can be characterized considering separately Case 1 and Case 2 (as introduced in Result 3.3).

Assuming Case 1:

For each n ∈ N , we define ρ s n as the value of the parameter ρ such that max

r * u ∈Ru r * u,n = r u n , ( 44
)
where r u n is defined by Eq. (29b) and the minimum value ρ s as

ρ s = min n∈N ρ s n . ( 45
)
The mitigation limit, denoted ρ ml , is defined as follows

ρ ml = ρ s , if ρ s < ρ lim (46a) ρ lim , if ρ s > ρ lim . (46b)

Assuming Case 2:

Because S has no fold points, ρ s n does not exist and therefore

ρ ml = ρ lim . ( 47
)
Results 3.3 and 3.4 are particularly interesting for engineering applications because they allow to discriminate harmful situations from harmless situations.

Analytical predictions proposed in Results 3.3 and 3.4 are illustrated and compared, for validation purposes, to numerical simulations in Sect. 5. The particular case of a linearized primary system coupled to N identical NESs is studied beforehand in next section.

Asymptotic analysis in the particular

case of a linear primary system coupled to N identical NESs

Literal expression of the mitigation limit

In this section, we assume N identical NESs (a n = a,

µ n = µ and α n = α, ∀n ∈ [1, N ]
) and a linear primary system ( fNL = 0). In this case, the following result can be obtained. ∀n ∈ [1, N ], max

r * ∈R r * n = w, ( 48
)
with w = 2 aρ + µρ (a 2 N -µρ) 3ρα . ( 49
)
The point defined as the constant coordinate w defines an unstable fixed point in the stable part of S.

Result 4.1 is proved in Appendix B. Eq. ( 49) shows that to ensure the existence of, at least, one nontrivial unstable fixed point of the slowflow (see Result 3.3), the expression (49) must be real, i.e. a 2 Nµρ > 0. We can therefore obtain the literal expression of ρ lim

ρ lim = a 2 N µ . ( 50
)
Moreover, because the NESs are identical, one also has

r u 1 = r u 2 = • • • = r u N -1 = r u N = r u = 2 √ 2 3 a + a 2 -3µ 2 α . (51)
Therefore, the limit value ρ s , solution of Eq. ( 44) for N identical NESs, can be obtained as

ρ s = a 2 µN 4a a 2 -3µ 2 + 5a 2 -3µ 2 (a 2 + µ 2 ) 2 . ( 52
)
One can be shown that the inequality ρ s < ρ lim is equivalent to µ < a/2. Therefore, according to Result 3.4, the literal expression of the mitigation limit for a linear primary system coupled to N identical NESs, denoted ρ * ml , is obtained and stated in the following result.

Result 4. [START_REF] Benoît | Chasse au canard ("duck hunting[END_REF] The mitigation limit ρ * ml , in the case of N identical NESs, is defined by the following equation

ρ * ml =            a 2 µN 4a a 2 -3µ 2 + 5a 2 -3µ 2 (a 2 + µ 2 ) 2 , if µ < a/2 (53a) a 2 N µ , if µ > a/2. ( 53b
)
Note that because a/2 < a/ √ 3, Eq. (53) is valid for both Case 1 and Case 2 introduced in Result 3.3.

Parametric study

In the case of identical NESs, the explicit expression of the mitigation limit (53) allows to easily perform a parametric study in order to obtain an optimal configuration, i.e. with the larger value of ρ ml . It is worth noting, regarding Eq. ( 53), that the mitigation limit is proportional to the number of NESs it only depends on the linear parameters of the NESs, namely a and µ.

Eq. ( 53) is plotted as a function of a and µ in Fig. 5, as a function of a for a fixed µ in Fig. 6(a) and as a function of µ for a fixed value of a in Fig. 6(b). For a given value of µ, ρ ml is an increasing function of a with the following asymptotic limit when a tends towards infinity lim a→+∞

ρ ml = 9N µ. ( 54
)
For a given value of a, it exists an optimal value of µ (denoted µ opt , the abscissas of the maximums in Figure 6(b)). The expression of µ opt can be obtained solving dρ ml /dµ = 0, leading to

µ opt = a 2 √ 3 -1. ( 55
)
Note that µ opt is always smaller than a/2. The corresponding mitigation limit is then computed substituting Eq. ( 55) into Eq. (53a) giving

ρ ml (µ opt ) = aN 2 9 + 6 √ 3. ( 56
)
It can be noted that the ρ ml (µ opt ) is proportional to aN . It might also be valuable to know, for a given mass of the mitigation system, if it is more effective to use one heavy NES or several light NESs. To answer the question Equation (56) (corresponding to N NESs of mass a) must be compared to the same substituting N by 1 and a by aN (corresponding to one NES of mass aN ) giving the same mitigation limit. Following this criterium, one heavy NES is equivalent to several light NESs.

Application to a Van der Pol oscillator

This section is devoted to illustrate and validate analytical predictions proposed in Sects. 3 and 4. To this end, a common Van der Pol oscillator is used as primary system. It is described by the following differential equation

ẍ + ρ ẋ x 2 -1 + x = 0. ( 57
)
corresponding to the following nonlinearity of the primary system f N L (x, ẋ) = ρx ẋ (see Eq. ( 3)). ρ is used as the bifurcation parameter. Since 0 < 1 and ρ > 0, the Van der Pol oscillator, described by Eq. ( 57), undergoes also relaxation oscillations and it can be shown that the LCO amplitude is A wo LCO = 2 (see for instance Example 5.2.1 in [START_REF] Kuehn | Multiple Time Scale Dynamics[END_REF]).

Performing the mathematical developments of Sect. 3, system of Eqs (37) reads

ρ H 1 (r 1 ) - H 2 1 (r 1 ) 4 - N k=1 r 2 k µ = 0 (58a) H n (r n ) -H 1 (r 1 ) = 0, for n = 2, . . . , N (58b)
Solving system of Eqs (58), fixed points of the slow-flow can be computed. Using Eq. ( 40) and Result 3.2, the stability of each fixed point can be determined.

Illustration of the results in the case of 3 different NESs

Let us first consider a case of a Van der Pol oscillator coupled to N = 3 different NESs using the following set of parameters

a 1 = 1 a 2 = 2 a 3 = 3 (59a) α 1 = 2 α 2 = 3 α 3 = 4 (59b) µ 1 = 0.3 µ 2 = 0.45 µ 3 = 0.6 (59c) = 0.001 (59d)
All NESs satisfy the condition µ n < a n / √ 3, therefore relaxation oscillations are possible in all planes (r n , s)

for n ∈ [1, 3].
The comparison between the theoretical bifurcation diagram, obtained from system of Eqs (58) and Result 3.2 and the maximum steady-state amplitudes obtained from numerical simulations of the non-averaged system [START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in Van der Pol-Duffing oscillator with nonlinear energy sink[END_REF] and of the slow-flow [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF], both including the Van der Pol primary nonlinearity, is presented in Figs. 7 and8.

The graphs of the maximum steady-state amplitude give a numerical estimation of the mitigation limit as the value of ρ for which the last jump of the amplitude is observed. In practice, we take the first value after the jump (with the larger amplitude). On Figs. 7 and8 this jump appears at ρ ≈ 9.3 for both the nonaverage system and the slow-flow. This value must be compared to the theoretical prediction ρ ml . The set of parameters used leads to a situation for which ρ ml is defined by Eq. (46a). Therefore, on Fig. 7, the mitigation limit corresponds to the smaller intersection of the branch of the larger fixed points and r u n . Here the smaller intersection is obtained for the third NES: we can read ρ ml ≈ 10 on Fig. 7(c), i.e. an overestimation of 6%. This is due to the fact that the prediction is obtained assuming 0 < 1, for example using = 10 -4 (resp. = 10 -2 ) we find an overestimation of 2% (resp. 19%), again for both the non-averaged system and the slow-flow.

Benchmark of the theoretical results in the case of identical NESs

In this section, the relevance of the obtained theoretical results is investigated in a systematic way. For this purpose, the following quantities are defined:

1. ρ n,1 ml : the first numerical estimation of the mitigation limit. It is "measured" on the graphs of the maximum steady-state amplitude obtained from numerical simulations of the non-averaged system [START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in Van der Pol-Duffing oscillator with nonlinear energy sink[END_REF]. ρ n,1 ml is the reference value of the mitigation limit.

2. ρ n,2 ml : the second numerical estimation of the mitigation limit. It is "measured" on the graphs of the maximum steady-state amplitude obtained from numerical simulations of the slow-flow [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF].

The previous numerical estimations of the mitigation limit are compared to the corresponding theoretical values:

1. ρ ml : the theoretical value of the mitigation limit defined in Sect. 3.4, considering a nonlinear primary system and valid for the general case of N different NESs.

2. ρ * ml : the theoretical value of the mitigation limit defined in Sect. 4.1, considering a linearized primary system. In this case a literal expression of the mitigation limit is obtained (see Eq. ( 53)) but it is valid only for N identical NESs.

The following parameters are used a n = 2 (60a)

α n = 4, ∀n ∈ [1, N ] (60b)
and the influence of the parameters µ n = µ (∀n ∈ [1, N ]), and N on the mitigation limit is investigated. Results are plotted in Figs. 9, 10 and 11 respectively. Relative errors of the mitigation limit are computed (in percentage) through the following expression

X -ρ n,1 ml ρ n,1 ml × 100, (61)
where X can denote ρ n,2 ml , ρ ml or ρ * ml . In the sequel, positive relative errors are called overestimations and negative relative errors are called underestimations.

In Fig. 9, one can see first that the slow-flow is able to estimate the mitigation limit. Indeed, the curves of ρ n,1 ml and ρ n,2 ml are almost superimposed with a maximum overestimation of 2.4% and 1.9% for N = 2 and N = 4 respectively. Then, the figure shows that the theoretical prediction ρ ml overestimates the mitigation limit when µ is close to its optimal value µ opt = 0.79 with a maximum overestimation of 11% and 18% for N = 2 and N = 4 respectively. However, the prediction is correct for µ > a/2, when the theoretical prediction is given by ρ s (see Result 3.4). Finally, we can observe that the theoretical prediction ρ * ml overestimates the mitigation limit for µ < a/2 (with a maximum overestimation of 7.4% and 14.3% for N = 2 and N = 4 respectively) and underestimates for µ > a/2 (with a maximum underestimation of 18.4% and 15.8% for N = 2 and N = 4 respectively). In Sect. 5.1, the differences between the reference value of the mitigation limit obtained with numerical simulations and the theoretical predictions have been explained by the assumption 0 < 1. Results presented in Fig. 9 show that neglecting the primary nonlinearity (to obtain ρ * ml ) is also a source of error, with a same order of magnitude.

In Fig. 10, the error made assuming 0 < 1 is investigated more precisely. We can see that the relative errors (61) of the theoretical predictions of the mitigation limit increase with the parameter . Table 1 presents the computed maximum and minimum relative errors. As expected, the maximum errors are obtained for = 10 -2 and the minimum errors for = 10 -4 , whatever the value of N . Finally, Fig. 11 shows the influence of the number N of parallel NESs. One can see that the overestimation of the theoretical predictions ρ ml and ρ * ml increase with N . We have minimum overestimations of 4.7% and 2.3% for ρ ml and ρ * ml for N = 1 respectively. The maximum overestimations are obtained for N = 6, reaching 18.7% and 14.9% for ρ ml and ρ * ml respectively. Observing Figs. 10 and11 and Table 1, it can be shown that the errors of the prediction ρ ml is often larger than the error of ρ * ml . This may seem contradictory, but in fact the error caused by the assumption of a small may be compensated by the error made ignoring the primary nonlinearity.

On the influence of the cubic stiffnesses of the NESs

The purpose of this final section is first to validate the separation method between harmless and harmful situations based only on the nature of the steady-state regimes. Then, in this context, the influence of the parameter α is investigated.

The maximum steady-state amplitudes obtained from numerical simulations of the slow-flow (10) including the Van der Pol primary nonlinearity as functions of µ and ρ are shown as density plots in Fig. 12 and as 3-D plots in Fig. 13 respectively for = 0.001, a = 2, N = 2 and three different values of the parameter α (2, 4 and 6). In Fig. 12, the mitigation limit ρ * ml (Eq. ( 53)) corresponding to α = 0 (linear primary system) is also drawn. In accordance with the order of magnitude of the nonlinearity of the primary system (see Eq. (3a)), the separation lines between small and large amplitudes for the three values of α show a good agreement with ρ * ml specially for small values of µ(< a/2). We observe that the amplitude jump corresponds actually to a large variation. The variation increases with α overall and the maximum values are obtained for small values of µ. For a given α the variation decreases when µ increases. For µ < a/2 the jumps occur between SMR and stable PR whereas for µ > a/ √ 3 they occur between stable PRs. The final but important remark is the fact that for large µ values and small α values there are no amplitude jumps anymore. This corresponds to situations in which the primary nonlinearity cannot be neglected compared to the nonlinearity of the NESs. It can be shown that the resulting bifurcation diagram of the slow-flow is composed of a single branch of stable fixed points.

Conclusion

We have performed an asymptotic analysis of a primary nonlinear unstable system coupled to N parallel Nonlinear Energy Sinks (NESs).

As usual, in the framework of NES properties exploration and particularly in the context of dynamic instabilities mitigation, four steady-state response regimes have been observed. They are classified into two cate-gories depending on whether the NESs act or not and therefore separating harmless situations from harmful situations. Three responses are referred to as harmless situations, namely complete suppression, partial suppression through periodic response, partial suppression through strongly modulated response; and one is referred to as harmful situation, that is called no suppression of the instability.

To analyze the steady-state response regimes, the system has been partitioned into slow-fast dynamics (denoted super-slow/slow dynamics for convenience) using the complexification-averaging method. The presence of a small dimensionless parameter related to the mass of the NES in the slow-flow system implies that it involves one "slow" complex variable and N "super-slow" complex variables. The "super-slow/slow" nature of the system allowed us to use an asymptotic approach to analyze it. In particular, the critical manifold of the slow-flow has been determined, it is a one dimensional parametric curve evolving in a N dimensional space. The shape of the critical manifold and the associated stability properties provide an analytical tool to predict the nature of the possible response Table 1: Relative errors of the theoretical predictions of the mitigation limit ρ ml and ρ * ml with respect to the reference ρ n,1 ml obtained using numerical simulations.

Minimum errors ( = 10 -4 ) Maximum errors ( = 10 ml (Eq. ( 7)) ρ n,2 ml (Eq. ( 10)) ml (Eq. ( 7)) ρ n,2 ml (Eq. ( 10)) and ρ n,2 ml and theoretical predictions ρ ml and ρ * ml of the mitigation limit as a function of µ using the set of parameters (60), with = 0.001 and for (a) N = 2 and (b) N = 4. regimes described above. In particular, the mitigation limit of the NESs, defined as the value of a chosen bifurcation parameter which separates harmful situations from harmless situations, has been predicted. Using more restrictive assumptions, i.e. neglecting the nonlinearity of the primary system and assuming N iden- ml (Eq. ( 7)) ρ n,2 ml (Eq. ( 10)) ml (Eq. ( 7)) ρ n,2 ml (Eq. ( 10))

(b)
(b)

Figure 10: Comparison between numerical estimations ρ n,1 ml and ρ n,2 ml and theoretical predictions ρ ml and ρ * ml of the mitigation limit as a function of using the set of parameters (60), with µ = µopt = 0.79 and for (a) N = 2 and (b) N = 4.

tical NESs, a literal expression has even been obtained. It has been highlighted that, in this context, the mitigation limit is proportional to the number of NESs and secondly that it depends only on the linear parameters of the NESs.

Finally, using a Van der Pol oscillator as a primary ml (Eq. ( 7)) ρ n,2 ml (Eq. ( 10))

Figure 11: Comparison between numerical estimation ρ n,1 ml and ρ n,2 ml and theoretical predictions ρ ml and ρ * ml of the mitigation limit as a function of N using the set of parameters (60) with µ = µopt = 0.79 and = 0.001. system, theoretical results have been compared, for validation purposes, to the numerical integration of the system. The comparison shows a good agreement as long as we remain within the limits of use of the asymptotic approach.

Compliance with Ethical Standards

Figure 1 :

 1 Figure 1: (a) Unstable linear dynamical system coupled to N parallel NESs. m is the mass, k is the stiffness and -p is the negative damping which causes the instability. (b) Zoom on the NESn.

Figure 2 :

 2 Figure 2: Example of typical critical manifold of a system containing two parallel NESs (N = 2) (a) in the (s, r1, r2)space (i.e. R + 3 ) and (b) projected on the (r1, r2)-plane (i.e. R + 2 ) . Following parameters are used: a1 = a2 = a3 = 1, α1 = 5, α2 = 6, α3 = 7, µ1 = 0.2, µ2 = 0.3 and µ3 = 0.4.

Figure 3 :

 3 Figure 3: Part of the Critical Manifold (CM) corresponding to the n-th NES, Eq (20). Following parameters are used: an = 1, αn = 6 and µn = 0.3.

Figure 4 :

 4 Figure 4: Examples of typical CM for a system containing (a) two and (b) three parallel NESs. Following parameters are used: a1 = a2 = a3 = 1, α1 = 5, α2 = 6, α3 = 7, µ1 = 0.2, µ2 = 0.3 and µ3 = 0.4.

Result 4 . 1

 41 Over the set of the fixed points R, the largest value for each n coordinate, with n ∈ [1, N ], is equal to w, i.e.

Figure 5 :

 5 Figure 5: Plot of Eq. (53) as a function of a and µ for N = 3.

Figure 6 :

 6 Figure 6: Plot of Eq. (53) as a function of (a) a with µ = 0.5 and (b) µ with a = 2 for N = 1, 2, 3, and 4.

Figure 7 :

 7 Figure 7: Comparison between the theoretical bifurcation diagram, obtained from system of Eqs (58) and Result 3.2 (magenta dots for stable fixed points and blue dots for unstable fixed points) and the maximum steady-state amplitudes obtained from numerical simulations of system of Eqs (7) (black empty circles) and system of Eqs (11) (green filled squares) including the Van der Pol primary nonlinearity. The bifurcation diagram and the maximum steady-state amplitudes are plotted for the variables (a) r1, (b) r2 and (c) r3. The set of parameters (59) is used.

Figure 8 :

 8 Figure 8: Bifurcation diagram and maximum steady-state amplitudes plotted for the variable s. Fig. (b) is a zoom of Fig. (a). The horizontal black line at s = 2 corresponds the LCO amplitude of the Van der Pol oscillator without any NES. The set of parameters (59) is used.

µ

  Mitigation limit (w.r.t. ρ)ρ * ml ρ ml ρ n,1

µ

  Mitigation limit (w.r.t. ρ)ρ * ml ρ ml ρ n,1

Figure 9 :

 9 Figure 9: Comparison between numerical estimations ρ n,1 ml

ǫ

  Mitigation limit (w.r.t. ρ)ρ * ml ρ ml ρ n,1

ǫ

  Mitigation limit (w.r.t. ρ)ρ * ml ρ ml ρ n,1

Figure 12 :Figure 13 :

 1213 Figure 12: Maximum steady-state amplitudes of the variable s obtained from numerical simulations of the slow-flow (10) including the Van der Pol primary nonlinearity as a function of µ and ρ for = 0.001, a = 2, N = 2 and (a) α = 2, (b) α = 4 and (c) α = 6, as density plots. The solid black lines with triangle markers represent the numerical estimation ρ n,2 ml

r 2 n

 2 = f 3 (r 2 1 ) = -r 2

  -2 )

		ρ ml	ρ * ml	ρ ml	ρ * ml
		N = 2 2.2% N = 4 4.6%	-1% 1.2%	31.9% 43.3%	27.7% 38.6%
	0	0.2 0.4 0.6 0.8

Conflict of interest

The authors declare that they have no conflict of interest concerning the publication of this manuscript.

A Projection of the super-slow dynamics on the critical manifold

Eq. ( 32) is developed into the following form ṡ = F (s, r 1 , . . . , r N , ϑ 1 . . . , ϑ N ) (62a) ṙn = G n (s, r 1 , . . . , r N , ϑ 1 . . . , ϑ N , ) , for n = 1, . . . , N (62b) θn = H n (s, r 1 , . . . , r N , ϑ 1 . . . , ϑ N , ) , for n = 1, . . . , N, (

where functions F, G n and H n are deduced from functions f and g n in Systems of Eqs. ( 10) and [START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF] as

G n = Re g n e -jθn (63b)

From system of Eqs (62) we obtain the real form of the slow and super-slow subsystems ( 13) and ( 14) as

and

We assume that fixed points of (62) exist only at the super-slow time scale and they are therefore fixed points of (65). Of course the CM can be also obtained solving Eqs. (65b) and (65c). Then, substituting system of Eqs [START_REF] Li | Dynamics of two vibro-impact nonlinear energy sinks in parallel under periodic and transient excitations[END_REF] with n = 1 (any n ∈ [1, N ] can be chosen) into (65a) , we obtain

which can be reduced to

where f r1 = 2 H 1 (r 1 )F. In reality, f r1 (r 1 , . . . , r N ) is a single-valued function because all variables r n (with

) are linked to each other through Eq. ( 22) (any n ∈ [1, N ] can be chosen as master component, we choose r 1 ).

B Proof of Result 4.1

Because N identical NESs are assumed, the function

Considering also a linear primary system ( fNL = 0), system of Eqs (37) becomes

Eq. (68b) is a third order polynomial equation with respect to r 2 n having the three following roots as a function of r 1

is the only real solution. Using system of Eqs (69), Eq. (68a) is written as follows

where g n (r 2 1 ) can be, for each term of the sum,

, Eq. (70) becomes

Ignoring the trivial solution, Eq. ( 71) is reduced to a second order polynomial equation with respect to r 2 1 with the following solutions

and

Considering first the case µ > a/ √ 3, r 2 n = f 1 (r 2 1 ) is the only real solution, r * 2 1,1 and r * 2 1,2 are therefore the only solutions given Eq. ( 48) with w = r * 1,2 .

We consider now the case for which µ < a/ √ 3. Solv-

1 ) (to find the intersection between the functions f 1 and f 2 ) and df 2 (x)/dx = 0 (to find the maximum of f 2 ) gives the same value of r 2 1 as

which corresponds also to one of the solutions of H (r 1 ) = 0 (see Eq. ( 25b)), i.e. r 1 = r m (where r m is r m n in the case of identical NESs). The difference between r * 2 1,2 and r 2 1 is

) and the right-hand side of Eq. ( 70) is bounded as follows

Let R l , R c and R r be the sets of fixed points r * , the solutions of H(r 1 ) = R l , H(r 1 ) = R c and H(r 1 ) = R r respectively. Consequently, because for r 2 1 > r 2 1 , H(r 1 ) in an increasing function and because of Eq. ( 76), the following inequalities hold

giving also Eq. ( 48) with w = r * 1,2 . To finish the proof, one must show that r * = [r 1 = w, . . . , r N = w] is an unstable fixed point on a stable part of S. To this end, the derivative with respect to r 1 of the function

defined by Eq. (39), is computed as