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

Scaling law for the slow flow of an unstable mechanical system coupled to a nonlinear energy sink

In this paper one first shows that the slow flow of a mechanical system with one unstable mode coupled to a Nonlinear Energy Sink (NES) can be reduced, in the neighborhood of a fold point of its critical manifold, to a normal form of the dynamic saddle-node bifurcation. This allows us to then obtain a scaling law for the slow flow dynamics and to improve the accuracy of the theoretical prediction of the mitigation limit of the NES previously obtained as part of a zeroth-order approximation. For that purpose, the governing equations of the coupled system are first simplified using a reduced-order model for the primary structure by keeping only its unstable modal coordinates. The slow flow is then derived by means of the complexificationaveraging method and, by the presence of a small perturbation parameter related to the mass ratio between the NES and the primary structure, it appears as a fast-slow system. The center manifold theorem is finally used to obtain the reduced form of the slow flow which is solved analytically leading to the scaling law. The latter reveals a nontrivial dependence with respect to the small perturbation parameter of the slow flow dynamics near the fold point, involving the fractional exponents 1/3 and 2/3. Finally, a new theoretical prediction of the mitigation limit is deduced from the scaling law. In the end, the proposed methodology is exemplified and validated numerically using an aeroelastic aircraft wing model coupled to one NES.

Introduction

The nonlinear energy sinks (NESs) are nowadays well-known devices used for passive shock and vibration mitigation of undesired oscillations caused by either external, parametric or self-excitations of a primary structure. An NES is classically defined as a nonlinear dynamical attachment consisting of a light mass (compared to the total mass of the primary structure), an essentially nonlinear spring (most of the time purely cubic) and a viscous linear damper. The operating of the NES is based on the phenomenon of targeted energy transfer (TET) through which a properly designed strongly nonlinear oscillator can be tuned to any frequency in order to perform an irreversible energy transfer from the primary structure towards itself. In their seminal papers [START_REF] Gendelman | Energy Pumping in Nonlinear Mechanical Oscillators: Part I-Dynamics of the Underlying Hamiltonian Systems[END_REF][START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators: Part II -Resonance capture[END_REF] Gendelman, Vakakis and coworkers explain the TET phenomenon by the interaction between two nonlinear modes of vibrations of the system producing a 1:1 resonance capture. Reviews of these concepts can be found in [3] and more recently in a part of [START_REF] Lu | Nonlinear dissipative devices in structural vibration control: A review[END_REF].

The use of NES to mitigate limit cycle oscillations (LCOs) stemming from dynamic instabilities has been widely studied in the past. The first work reported in this framework concerns the mitigation of LCOs of the Van der Pol oscillator [START_REF] Lee | Suppression of limit cycle oscillations in the van der Pol oscillator by means of passive non-linear energy sinks[END_REF]. By means of a perturbation analysis, Gendelman and Bar [START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] predicted the response regimes of the same system. This work has been extended to a Van der Pol-Duffing oscillator coupled to one NES by Domany and Gendelman [START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in van der pol-duffing oscillator with nonlinear energy sink[END_REF]. Numerous papers have been dedicated to the problem of mitigation, by means of one or several NESs, of LCOs due to aeroelastic instabilities. The seminal papers are those of Lee et al. [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] Lee | Enhancing the robustness of aeroelastic instability suppression using multi-degree-of-freedom nonlinear energy sinks[END_REF] and concerned aeroelastic instabilities in aircraft wing. In these works, the problem has been investigated both numerically and experimentally. The same aircraft wing model coupled to one NES has been studied theoretically using the so-called complexification-averaging method [START_REF] Manevitch | Complex representation of dynamics of coupled nonlinear oscillators[END_REF] together with a singular perturbation approach [START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF] or a perturbation algorithm based on a mixed multiple scale-harmonic balance method (MSHBM) [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF]. Vaurigaud et al. [START_REF] Vaurigaud | Passive control of aeroelastic instability in a long span bridge model prone to couple dflutter using targeted energy transfer[END_REF] investigated a problem of passive nonlinear TET between a two degrees of freedom long span bridge model prone to coupled flutter and an NES. The mitigation, using an NES, of vortex-induced vibrations resulting from the nonlinear interaction of a laminar flow and a rigid circular cylinder has been first studied by Tumkur et al. [START_REF] Tumkur | Reducedorder model for laminar vortex-induced vibration of a rigid circular cylinder with an internal nonlinear absorber[END_REF] constructing a two-DOF reduced-order model of the system and validating it by means of a comparison with a full order finite-element model. On the same issue, an improved and experimentally validated reduced model has been proposed by Dai et al. [START_REF] Dai | Vortex-induced vibrations mitigation through a nonlinear energy sink[END_REF]. The passive control of helicopter ground resonance instabilities by a means of NESs has been theoretically and numerically studied by Bergeot et al. considering an NES attached to the helicopter fuselage [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF] and to the blades [START_REF] Bergeot | Passive suppression of helicopter ground resonance using nonlinear energy sinks attached on the helicopter blades[END_REF]. Chatter control in machine tool vibrations has been studied by Gourc et al. [START_REF] Gourc | Chatter Control in Turning Process with a Nonlinear Energy Sink[END_REF] considering a vibro-impact NES and by Nankali et al. [START_REF] Nankali | Targeted energy transfers for suppressing regenerative machine tool vibrations[END_REF] considering a purely cubic NES. 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], mitigation of friction-induced vibrations due to mode coupling instabilities in a friction system has been investigated. The use of several NESs to mitigate LCOs has been analytically investigated by Bergeot et Bellizzi. First, a network of parallel NESs coupled to a Van der Pol oscillator has been analyzed in [START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF] and then the prediction of the steady-state regimes of a multi-DOF dynamical system having one unstable mode and coupled to a set of NESs has been performed in [START_REF] Bergeot | Steady-state regimes prediction of a multi-degree-of-freedom unstable dynamical system coupled to a set of nonlinear energy sinks[END_REF]. Finally, in a recent paper, the possibility of mitigating simultaneously two unstable modes of a linear multi-DOF primary system using a single NES has been investigated by Bergeot et al. [START_REF] Bergeot | Dynamic behavior analysis of a mechanical system with two unstable modes coupled to a single nonlinear energy sink[END_REF]. The asymptotic analysis proposed in the latter reference reveals that the NES attachment can produce several and complex mitigated responses which results from the presence of several stable solutions and from the competition between them.

In general, when an NES is attached to a primary structure, the resulting coupled model is analyzed by introducing a small perturbation parameter related to the mass ratio between the NES and the primary structure. Under the assumption of a 1:1 resonance capture, the system is first averaged over one period corresponding to a natural frequency of the primary structure. The resulting averaged system, called slow flow, is then analyzed by means of singular perturbation techniques (multiple scales method [START_REF] Nayfeh | Perturbation Methods[END_REF] or geometric singular perturbation theory [START_REF] Jones | Geometric singular perturbation theory[END_REF]). The first key point of these analytical treatments is the partition of the slow flow into two time scales (one slow and one fast 1 ). In this representation the time evolution of the slow flow is thus described as a succession of slow and fast epochs which are described independently. The second key point is the definition of the critical manifold whose system trajectories approaches during slow epochs. In general, the analysis is carried out within the zeroth-order approximation, i.e. within the limit case in which the perturbation parameter is equal to zero. Consequently, the obtained analytical results depreciate for the largest values of the perturbation parameter by thus limiting their predictive power. Luongo and Zulli [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF] used a perturbation technique to express the invariant manifold as a power series with respect to the perturbation parameter in which the first term is the critical manifold2 . As stated by the Fenichel theory [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF], this representation is relevant for normally hyperbolic branch of the critical manifold but fails, at fold points (the junction between slow and fast epochs takes place in the vicinity of these points), to approach the actual trajectory of the system. Indeed, the literature on dynamical systems shows that the trajectory passes through the region near these fold points with a nontrivial scaling behavior with respect to the perturbation parameter that it is not possible to obtain with the classical perturbation methods (see e.g. [START_REF] Christian | Multiple Time Scale Dynamics[END_REF] Sect. 5.4).

By means of a sophisticated analysis, the originality of the present work is to provide an analytical description of the system dynamic behavior which reports on this nontrivial dependence on the small perturbation parameter. This analytical result is hereafter referred as the scaling law of the slow flow, this notion is clarified below within the main body of the present manuscript. The main objective is to obtain an analytical description of the system with a high predictive power compared to what has been proposed in the past in the literature and then to hope in the future of including them in a NES design procedure.

The equations of motion of the system under study are presented in Section 2. The full order model equations of a mechanical system coupled to one NES are derived in Section 2.1. Assuming a primary structure with only one unstable mode, the modal reduction technique, leading to the reduced-order dimensionless model, is presented in Section 2.2. Section 3 summarizes and extends some previous analytical results of the literature obtained within the zeroth-order approximation. Section 4 expounds the main results of the paper. In Section 4.1, the center manifold theorem is used to reduce the dynamics of the slow flow of the system. This reduced form is solved analytically in Section 4.2 leading to the scaling law from which a new analytical prediction of the mitigation limit is derived in Section 4.3. In Section 5, the proposed methodology is exemplified and validated numerically using an aeroelastic aircraft wing model coupled to one NES. Finally, concluding remarks are formulated in Section 6.

Equations of the model

The dynamical system under study in this paper is presented in the present section. The following developments are similar to those presented in [START_REF] Bergeot | Steady-state regimes prediction of a multi-degree-of-freedom unstable dynamical system coupled to a set of nonlinear energy sinks[END_REF].

The initial model

The primary structure under consideration in this work is described by the following general set of nonlinear differential equations

Mẍ + C ẋ + Kx + gNL (x) = 0, (1)
where x = (x 1 , . . . , xN ) T with () T the transpose operator, the dot represents time-differentiation and M, C and K are the mass matrix, the damping matrix and the stiffness matrix respectively. The nonlinear vector-valued function gNL represents the nonlinearity of the primary system. Moreover, cubic nonlinearities are assumed, i.e. each component of gNL is a linear combination of monomial terms of order 3. In these first equations, the tilde symbol is used before the rescaling which will be carried out below (to obtain variables without tilde after rescaling). System (1) can undergo a dynamic instability of the trivial solution through a Hopf bifurcation, i.e. it has a pair of complex conjugate eigenvalues which cross the complex plane imaginary axis. The loss of stability of the trivial solution comes with the production of a periodic solution and the nonlinear function gNL allows the existence of Limit Cycle Oscillations (LCOs) on which the system can saturate. For the purpose of mitigating the LCOs, one purely cubic ungrounded NES with mass mh , damping coefficient ch and cubic stiffness kNL h is used. Taking into account the NES displacement h(t), the equations of the coupled system are

Mẍ + C ẋ + Kx + gNL (x) + B ch A ẋ -ḣ + kNL h Ax -h 3 = 0 (2a) mh ḧ -ch A ẋ -ḣ -kNL h Ax -h 3 = 0, (2b)
where A = (a 1 , . . . , a N ) and B = b1 , . . . , bN T are the influence coefficient vectors which depend on the position of the NES.

One assumes now that the mass mh and the damping coefficient ch of the NES are small compared to the mass of the primary structure (1) and also that the latter is weakly nonlinear. We therefore introduce a small perturbation dimensionless parameter (with 0 < 1) and the associated rescaled coefficients m h and c h and vector function g NL as mh = m h , ch = c h and gNL = g NL . Then, the variables x and h are also rescaled through as x = x √ and h = h √ . Due to the latter rescaling and because cubic nonlinearities are assumed one has g NL (x) = 3/2 g NL (x). Inserting m h , c h , g NL , x and h into (2), multiplying Eq. (2a) by M-1 and Eq. (2b) by 1/( m h ) and finally neglecting terms of order strictly higher than 1 in lead to

ẍ + C ẋ + Kx + B ḧ = 0 (3a) ḧ -η h (A ẋ -ḣ) -k NL h (Ax -h) 3 = 0. ( 3b
)
with

C = M-1 C, K = M-1 K, B = m h M-1 B, η h = c h /m h and k NL h = kNL h /m h .
Note that because of the previous assumptions and rescaling, the nonlinearities of the primary structure are neglected from now. Therefore, nonlinear modal interactions in the primary structure such as internal resonances are not taken into account in this work.

Reduction of the dynamics

The system (3) is now reduced taking into account that the primary system (1) undergoes a single instability of the trivial solution. The objective is to obtain a system with two degrees of freedom (one for the unstable mode of the primary structure and one for the NES) on which one will be able to carry out an analysis which would be impossible on a the original system.

To achieve that, it is first convenient to introduce new coordinates v = x + Bh and w = h-Ax giving reciprocally

x = (I + BA) -1 (v -Bw) ≈ (I -BA)(v -Bw) (4a) h = (1 + AB) -1 (w + Av) ≈ (1 -AB)(w + Av) ( 4b
)
where I is the identity matrix of size N . Using Eq. ( 4), Eq. ( 3) is transformed into the following form

v + C v + Kv -CB (A v + ẇ) + KB (Av + w) = 0 (5a) ẅ + η h ẇ + k NL h w 3 -AC v -AKv+ ABη h ẇ + ACB (A v + ẇ) + AKB (Av + w) + ABk NL h w 3 = 0. ( 5b
)
where only terms of order lower or equal to 1 in were kept.

To retain only the essential features of the single instability, the reduction of the dynamics is performed on Eq. (5a). To achieve this, the latter, is written in state-space form as follows

ẏ = Dy + [D 1 y + D 1 w + D 3 ẇ] ( 6
)
where y = (v 1 , . . . , v N , v1 , . . . , vN ) T ,

D = 0 I N -K -C , D 1 = 0 0 KBA CBA , D 2 = 0 KB and D 3 = 0 CB . ( 7
)
Because the matrix D, which characterizes the linear dynamics of the primary system, is not symmetric, its diagonalization must be performed using biorthogonality property of the right eigenvectors r i (i = 1, . . . , 2N ) and the left eigenvectors l j (j = 1, . . . , 2N ) of D. When the eigenvectors are properly normalized, the biorthogonality properties consist in L T DR = Λ where

R = [r 1 r 1 . . . r N r N ], L = l 1 l 1 . . . l N l N and Λ = diag(λ 1 , λ 1 , . . . , λ N , λ N )
are the right and left eigenvector matrices and the eigenvalue diagonal matrix respectively.

The so-called biorthogonal coordinates are then introduced. They are constituted of N pairs of complex conjugates, q n and q n (n = 1, . . . , N ), and defined by the following relations y = Rq and q = L T y, [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF] where q = (q 1 , q 1 , . . . , q N , q n ) T . Substituting Eq. (8) into Eq. ( 6), the equations of motion take the form of the following system

q = Λq + L T [D 1 Rq + D 2 w + D 3 ẇ] , (9)
and therefore Eq. ( 5) is equivalent to

q -Λq -L T [D 1 Rq + D 2 w + D 3 ẇ] = 0 (10a) ẅ + η h ẇ + k NL h w 3 - AC(R DL q U + R DR q D ) -AK(R UL q U + R UR q D )+ ABη h ẇ + ACB A(R DL q U + R DR q D ) + ẇ + AKB A(R UL q U + R UR q D ) + w + ABk NL h w 3 = 0, ( 10b
)
where the matrix R and the vector q have been split into a N × N -block matrix and a N × 1-block respectively as

R = R UL R UR R DL R DR
and q = ((q U ) T , (q D ) T ) T [START_REF] Manevitch | Complex representation of dynamics of coupled nonlinear oscillators[END_REF] where "U", "D", "L" and "R" in previous superscripts must be understood as the initial letters of up, down, left and right respectively.

As previously mentioned, one mode of the primary system can become unstable through a Hopf bifurcation. Without loss of generality, the first mode is chosen to be the potential unstable mode. That means that when the chosen physical bifurcation parameter, denoted σ, crosses the particular parameter value σ Hopf , called Hopf bifurcation point, a pair of complex conjugate eigenvalues crosses the complex plane imaginary axis (without loss of generality this is first pair (λ 1 , λ 1 )) whereas all the other eigenvalues (λ n , λ n ) (with n = 2, . . . , N ) have strictly negative real parts.

Neglecting the coupling term L T [D 1 Rq + D 2 w + D 3 ẇ], Eq. (10a) reduces to the diagonalized system q -Λq = 0. Therefore, in this case, the variables q n and q n (n = 2, . . . , N ), related to the stable modes, tend to zero. When the coupling term (which is of the order of magnitude of ) is taken into account, the variables q n and q n do not vanish but their contributions are very small and can be neglected [START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF][START_REF] Bergeot | Steady-state regimes prediction of a multi-degree-of-freedom unstable dynamical system coupled to a set of nonlinear energy sinks[END_REF][START_REF] Bergeot | Dynamic behavior analysis of a mechanical system with two unstable modes coupled to a single nonlinear energy sink[END_REF]. Consequently, all terms related to q n and q n (with n = 2, . . . , N ) are omitted from further consideration and Eq. ( 9) is reduced to

q1 = λ 1 q 1 + l T 1 D 1 (r 1 q 1 + r 1 q 1 ) + D 2 w + D 3 ẇ . ( 12
)
Finally, grouping Eqs. (10b) and ( 12), splitting λ 1 into its real and imaginary parts as

λ 1 = λ R + jλ I (with j 2 = -1) and switching the time scale as t → λ I t yield q1 - λ R λ I + j q 1 - λ I l T 1 D 1 (r 1 q 1 + r 1 q 1 ) + D 2 w + D 3 ẇ = 0 (13a) ẅ + η h λ I ẇ + k N L h λ 2 I w 3 - 1 λ 2 I AC(r DL 1 q 1 + r DL 1 q 1 ) + AK(r UL 1 q 1 + r UL 1 q 1 ) + AB η h λ I ẇ + 1 λ 2 I ACB A(r DL 1 q 1 + r DL 1 q 1 ) + λ I ẇ + 1 λ 2 I AKB A(r UL 1 q 1 + r UL 1 q 1 ) + w + AB k N L h λ 2 I w 3 = 0. ( 13b
)
where the vectors r UL 1 and r DL 1 correspond to the first columns of the matrices R DL and R UL respectively (see Eq. ( 11)).

The system of equations ( 13) is the system under study in the following mathematical developments.

Asymptotic analysis: the zeroth-order approximation of the slow flow

The present section summarizes and extends previous results of the literature [START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF][START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF][START_REF] Bergeot | Steady-state regimes prediction of a multi-degree-of-freedom unstable dynamical system coupled to a set of nonlinear energy sinks[END_REF]. In particular, an analytical expression of the mitigation limit (hereafter called the zeroth-order approximation of the mitigation limit) is derived at the end of the section under several assumptions. The first is given just below.

In Eq. ( 13), the quantities λ R , λ I , the vectors l T 1 , r 1 , r UL 1 and r DL 1 and the matrices C, K (and those which depend on them) depend in general on the bifurcation parameter under consideration σ. In this work one assumes that the NES effects are not too far from the Hopf bifurcation point σ Hopf and therefore the latter are evaluated at σ Hopf (using the notation () Hopf ) expect λ R which is by definition equal to zero at σ Hopf . λ R is evaluated at the actual value of σ and the dimensionless eigenvalue real part λ R /ω (where ω = λ I,Hopf , the imaginary part of λ 1 evaluated at σ Hopf ) is supposed to be small (i.e. of the order of O( )). The parameter ρ is then finally introduced as

ρ = λ R ω ( 14
)
where ρ is of the order of O(1).

The rescaled dimensionless eigenvalue real part ρ of the unstable mode is considered as a generalized bifurcation parameter (that results from the physical bifurcation parameter σ of the initial mechanical system) from which one will follow the changes in the system dynamic behavior.

Governing equations of the slow flow

It is know [START_REF] Gendelman | Energy Pumping in Nonlinear Mechanical Oscillators: Part I-Dynamics of the Underlying Hamiltonian Systems[END_REF][START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators: Part II -Resonance capture[END_REF] that TET appears in the vicinity of the 1:1 resonance between two nonlinear modes of the coupled structure, the latter resulting from the interaction between a linear mode of the primary structure (here this is the considered unstable mode) and the nonlinear mode of the NES. This phenomenon is called 1:1 resonance capture and occurs at a frequency close to the natural frequency of the linear primary structure. To investigate the solution in the neighborhood of this 1:1 resonance capture the complexification-averaging method [START_REF] Manevitch | Complex representation of dynamics of coupled nonlinear oscillators[END_REF]3] is used. The resulting complex averaged dynamics is called slow flow.

The first step consists in introducing the following complex variables

ψ = ẇ + jw ( 15
)
or equivalently

w = ψ -ψ 2j , ẇ = ψ + ψ 2 and ẅ = ψ - j 2 ψ + ψ ( 16
)
and expressing the complex variables q 1 and ψ as

q 1 = φ1 e jt and ψ = φ 2 e jt . ( 17
)
Substituting ( 17) into [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF] and averaging over the period T = 2π leads to the following system of complex slow modulated amplitudes

φ1 = ρ + l T 1,Hopf D 1,Hopf r 1,Hopf φ1 + 1 2 -jl T 1,Hopf D 2,Hopf + l T 1,Hopf D 3,Hopf φ 2 (18a) φ2 = AK Hopf r UL 1,Hopf + AC Hopf r DL 1,Hopf ω 2 φ1 - 1 2 j + µ -j 3α 4 |φ 2 | 2 φ 2 - 1 ω 2 AK Hopf BTr UL 1,Hopf + 1 2 AC Hopf BTr DL 1,Hopf φ1 + 1 2 AC Hopf B ω -j AK Hopf B ω 2 + AB µ -j 3α 4 |φ 2 | 2 φ 2 (18b) with µ = η h ω and α = k N L h ω 2 (19)
and recalling that the notation () Hopf means that the considered quantity is evaluated at the Hopf bifurcation point σ Hopf .

Introducing the variable φ 1 = 2 AK Hopf r UL 1,Hopf + AC Hopf r DL 1,Hopf φ1 /ω 2 , and ignoring the terms of the order of magnitude of in (18b) 3 , Eq. ( 18) reduces to the following synthetic form

φ1 + [(a -ρ) φ 1 + bφ 2 ] = 0 (20a) φ2 + 1 2 φ 1 + 1 2 (µ + j)φ 2 + j 3 8 αφ 2 |φ 2 | 2 = 0. ( 20b
)
with the complex coefficients a and b defined as

a = -l T 1,Hopf D 1,Hopf r 1,Hopf b = 1 ω 2 -jl T 1,Hopf D 2,Hopf + l T 1,Hopf D 3,Hopf AK Hopf r UL 1,Hopf + ACr DL 1,Hopf . (21)
Then, substituting φ 1 and φ 2 in Eq. ( 18) by their polar coordinates defined as φ 1 = re jθ 1 and φ 2 = se jθ 2 , new equations of motion for the real amplitudes r and s and phase difference ∆ = θ 2 -θ 1 can be obtained as

ṙ = f (r, s, ∆) (22a) ṡ = g1 (r, s, ∆, ) (22b)
∆ = g2 (r, s, ∆, ) (22c)
with

f (r, s, ∆) = (ρ -a R ) r -(b R cos ∆ -b I sin ∆) s, (23a) g1 (r, s, ∆, ) = - 1 2 (µs + r cos ∆) , (23b) g2 (r, s, ∆, ) = 3 8 αs 2 - 1 2 + r 2s sin ∆ + a I + b I cos ∆ r + b R sin ∆ r (23c)
where the complex parameters a and b were split into real and imaginary parts as

a = a R + ja I and b = b R + jb I .
Because of the presence of the small parameter , the slow flow is governed by different time scales. More precisely, within the framework of the geometric singular perturbation theory [START_REF] Jones | Geometric singular perturbation theory[END_REF], Eq. ( 22) appears as a (2, 1)-fast-slow system where s and ∆ are the fast variables and r the slow variable. The latter is related to the modal component q 1 whereas the fast variables are associated with the relative displacement w between the NES and the primary structure.

Critical manifold and fold points of the slow flow

The real form of the slow flow, i.e. Eq. ( 22), is written with respect to the dimensionless slow time

τ = t as follows r = f (r, s, ∆) (24a) s = g1 (r, s, ∆, ) (24b)
∆ = g2 (r, s, ∆, ) (24c)
where (.) denotes the derivative with respect to the slow time τ . Considering = 0 yields the slow subsystem

r = f (r, s, ∆) (25a) 0 = g1 (r, s, ∆, 0) (25b) 0 = g2 (r, s, ∆, 0) (25c)
which is a differential-algebraic equation, and the fast subsystem

ṙ = 0 (26a) ṡ = g1 (r, s, ∆, 0) (26b) ∆ = g2 (r, s, ∆, 0) (26c)
The critical manifold of the slow flow is the solution of the algebraic part of Eq. ( 25) and it is expressed as follows

M 0 = (r, s, ∆) ∈ R + 2 × [-π, π] g1 (r, s, ∆, 0) = 0 and g2 (r, s, ∆, 0) = 0 . ( 27
)
Combining Eqs. (25b) and (25c) leads to the following modulus and argument equations

r = s µ 2 + 1 - 3αs 2 4 2 = H(s) (28a) tan ∆ = 3αs 2 -4 4µ . ( 28b
)
Stability analysis of the critical manifold is now carried out considering the fast subsystem [START_REF] Jones | Geometric singular perturbation theory[END_REF]. It can be shown that a fixed point of ( 26) is stable if d s H(s) > 0 and unstable if d s H(s) < 0. Hence, the subset of M 0 satisfying d s H(s) < 0 defines the attracting zone for the fast dynamics whereas the subset of M 0 satisfying d s H(s) > 0 define the repelling zone.

Exploiting the polynomial properties of H, it can be shown that the local extrema of H (i.e

d s H(s) = 0) occur at s LF = 2 3 √ α 2 -1 -3µ 2 and s RF = 2 3 √ α 2 + 1 -3µ 2 (29)
if the following relation holds

µ < 1 √ 3 . ( 30
)
In Eq. ( 29), s LF and s RF are the abscissa values in the (s, r)-plane of the a maximum and the minimum of the function H respectively. The superscripts () LF and () RF refer to left fold point and right fold point respectively. Indeed, in the (s, r, ∆)-space, the two points (s RF , r RF , ∆ RF ) and (s LF , r LF , ∆ LF ) (where r LF , r RF , ∆ LF and ∆ RF can be deduced from s LF and s RF using Eq. ( 28)) are generally called fold points.

When the condition [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF] is not satisfied the function H no longer has local extrema. In the remaining of the paper, one considers situations in which (30) holds.

A typical critical manifold (more precisely its real part defined by Eq. (28a)) is depicted in Fig. 1. One can see that the fold points connect attracting parts to repelling part of the critical manifold M 0 . The two scalars s D and s U , which are the horizontal projection of the fold points on the critical manifold, are defined by

H s RF = H s D and H s LF = H s U giving s D = 2 √ 2 3 √ α 1 -1 -3µ 2 and s U = 2 √ 2 3 √ α 1 + 1 -3µ 2 . ( 31
)

Fixed points and fold singularities of the slow flow

Because 0 < 1, the fixed points of the slow flow [START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF], which characterize periodic solutions of [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF], can be approximated by those of the slow subsystem [START_REF] Nayfeh | Perturbation Methods[END_REF]. From that zeroth-order approximation an analytical expression of these fixed points can be obtained. Note that within this approximation, the fixed points are necessarily located on the critical manifold.

First, substituting Eq. (28a) into Eq. (25a) results in (H(s)) = f (H(s), s, ∆) which reduces to

d s H(s)s = f s (s), (32)
where the explicit expression of the function f s is not given here. Then, a (regular) fixed point of Eq. ( 32) is defined as the roots of f s (s) = 0 with d s H(s) = 0. Solving f s (s) = 0 gives two fixed points (in addition to the trivial solution) whose analytical expressions are

s e 1 = 2 3α (a R -ρ) 2a R -b I -2ρ + b 2 I -4µ (a R -ρ) (µa R -b R -µρ) 1/2 , ( 33a
)
s e 2 = 2 3α (a R -ρ) 2a R -b I -2ρ -b 2 I -4µ (a R -ρ) (µa R -b R -µρ) 1/2 , ( 33b
)
with s e 1 < s e 2 . As usual, the stability of the fixed points are found looking for the sign of d s fs(s) H (s) . Both fixed points s e 1 and s e 2 no longer exist when they become complex for

ρ S = 2µa R + b 2 I + b 2 R -b R 2µ ( 34
)
where the superscript () S has no particular meaning here, it is just there to give a name to this particular value of ρ.

A singular fixed points of Eq. ( 32), also called folded singularites, are defined as the roots of the following nonlinear equations f s (s) = 0 and d s H(s) = 0. Folded singularities are hints of particular solutions of such fast-slow systems, called canard cycles [START_REF] Benoît | Chasse au canard ("duck hunting[END_REF]. These solutions are not investigated in this paper.

The possible responses

Previous analysis provides a qualitative description of the slow flow dynamics in the limit case where = 0, i.e. within a zeroth-order approximation. Within this approximation, the slow flow evolves at two time scales: a slow time scale in which the slow flow is on the critical manifold and described by the slow subsystem [START_REF] Nayfeh | Perturbation Methods[END_REF] and a fast time scale in which the slow flow is outside the critical manifold and described by the fast subsystem [START_REF] Jones | Geometric singular perturbation theory[END_REF]. The particular S-shape of the critical manifold together with the stability analysis of the fixed points of the slow flow allows us to explain and predict its different responses and consequently those of the initial full order system [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF]. As has been widely discussed in the literature (see e.g. [START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a 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]), four scenarios are possible. In previous works by the authors [START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF][START_REF] Bergeot | Steady-state regimes prediction of a multi-degree-of-freedom unstable dynamical system coupled to a set of nonlinear energy sinks[END_REF] these responses are classified into two categories. In the first category, the NES acts, resulting in the three responses called harmless situations. First, the Complete suppression in which the trivial fixed point of the slow flow ( 22) is stable and then reached. In this case, the system is stabilized because of the additional part of the NES including mass and damping. In general, because a light mass and a small damping are considered for the NES, this stabilization effect is negligible, in term of bifurcation parameter range, compared to the mitigation effects described hereinafter. In other words, the Hopf bifurcation point of the primary structure without NES is usually very slightly smaller than the Hopf bifurcation point of the coupled system. Of course, in a context of nonlinear vibrations absorption, the complete suppression is not the desired effect. Then, one can observe mitigation through periodic response.

In this case, a nontrivial stable fixed point of the slow flow ( 22) is reached, that corresponds to a periodic regime for the initial system [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF]. The last harmless situations correspond to mitigation through Strongly Modulated Responses (SMRs). SMRs correspond to relaxations oscillations of the slow flow [START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF], that corresponds to a quasi-periodic (amplitude and phase modulated) regime for the initial system [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF]. In the second category, called harmful situation, the NES is not able to produce small amplitude responses. For the initial nonlinear system one observes in this case a limit cycle with an amplitude close to that of the system without NES. Note that if the nonlinearity of the primary structure is neglected in the governing equations of the slow flow (this the case in the present work), the harmful situation corresponds for it to an unbounded regime in which its trajectory growths to infinity.

In this paper, we are particularly interested in the transition from SMR to no mitigation which is also, in general, the transition from harmless to harmful situations. Therefore one recalls here briefly what relaxation oscillations of the slow flow are. Note that the following description holds in the limit case for which = 0. In the (s, r)-plane, from an initial condition near zero and outside the critical manifold, the trajectory of the slow flow evolves rapidly and horizontally to the left attracting branch of M 0 (this fast epoch is described by the fast subsystem ( 26)). Then, the slow flow evolves slowly on this branch (this slow epoch is described by the slow subsystem ( 25)) and, if there is no stable fixed points on it, the trajectory reached the left fold point (s LF , r LF ) at which M 0 becomes repelling. From the fold point, the slow flow undergoes a fast horizontal jump to (s U , r LF ) on the right attracting part of M 0 and then a second slow epoch to the right fold point (s RF , r RF ). A second horizontal jump occurs from (s RF , r RF ) to (s D , r RF ) and the trajectory returns to the left attracting part of the M 0 . A third slow epoch occurs on M 0 to (s LF , r LF ) and so on.

Zeroth-order analytical prediction of the mitigation limit

The mitigation limit is defined below. Definition 3.1. Considering a set of initial conditions (for the slow flow) as a small perturbation of the trivial solution, the mitigation limit is defined as the value of the generalized bifurcation parameter ρ which separates harmless situations from harmful situation.

As discussed in the literature [START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a 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][START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF], the mitigation limit (denoted here ρ * 0 ) is defined, within the zeroth-order approximation, as the value of ρ for which the fixed point s e 2 becomes smaller than s U defined by Eq. ( 31) or for which s e 2 no longer exists (see Eq. ( 34)). Indeed, when the largest unstable fixed exists (with, in addition, s e 2 > s U ), it prevents the phase trajectory to grow to infinity on the right branch of the critical manifold.

From Eqs. ( 31) and (33a) the solution of s e 2 = s U is derived as

ρ U 0 = a R + 2 1 -3µ 2 b I + b I -2 1 -3µ 2 + 7 µ 2 b I (µ 2 + 1) 2 + 3µ 2 -4 1 -3µ 2 -5 µb R (µ 2 + 1) 2 . ( 35
)
Then using (34), the following conditional expression for ρ * 0 is obtained:

ρ * 0 = ρ U 0 (Eq. (35)), if µ < µ * (36a) ρ S (Eq. (34)), if µ > µ * , ( 36b
)
with µ * the special value of µ solution of ρ U 0 = ρ S which is in general smaller than 1/ √ 3.

We can notice that within the zeroth-order approximation, the mitigation limit does not depend on the nonlinear parameter α of the NES but only its damping coefficient µ (together with the parameters of the primary system). Therefore, for given parameters of the primary structure, it exists an optimal value of µ, denoted µ opt 0 , which maximizes the mitigation limit. Consequently, assuming that µ opt 0 < µ * , µ opt 0 is the solution of ∂ µ ρ U 0 = 0. Solving the latter equation yields

µ opt 0 = - 1 B b R b I + B |b| 4 b 2 R 1/3 -3b 2 R b 2 I + b 2 R --B |b| 4 b 2 R 1/3 -2Bb 2 R + 2b 2 I b 2 R - 2b I b R b 2 I + b 2 R 8b 2 I + 9b 2 R B |b| 4 b 2 R 1/3 -3b 2 R b 2 I + b 2 R 1/2 (37)
with B = 3b 2 R + 4b 2 I and where the subscript 0 is used to highlight that the expression (37) is obtained within the zeroth-order approximation.

Scaling law for the slow flow near a fold point

The Fenichel theory [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] guarantees the persistence of the relaxation oscillations scenario, described at the end of Section 3.4, for 0 < 1 by stating that Eq. ( 24) has an invariant manifold M in the O( )-vicinity of M 0 and with the same stability properties with respect to the fast variables as M 0 (attracting or repelling). In this case, concerning the first jump described previously, the jump point at which the trajectory leaves the attracting branch of M to undergo the jump is no longer the fold point (s LF , r LF ). The same is true for the arrival point at which the trajectory reaches the right attracting part of M after the jump which is no longer (s U , r LF ). The goal of this section is to determine these jump and arrival points and consequently to obtain an analytical prediction of the mitigation limit which, unlike Eq. ( 36), takes into account .

To achieve that, one must determine beforehand the scaling law of the slow flow dynamics in the neighborhood of the fold point (s LF , r LF ), i.e. a law which describes the -dependance of the distance between M 0 and the actual trajectory of the slow flow for 0 < 1. The jump and arrival points are then deduced from this scaling law.

The Fenichel theorem is true for normally hyperbolic branch of M 0 (i.e. outside the fold points) and fails to approach the actual trajectory of the system at fold points. Indeed, the literature on dynamical systems shows that the trajectory still passes through the region near these fold points but with a nontrivial scaling behavior involving the fractional exponents 1/3 and 2/3 for the perturbation parameter (see e.g. [START_REF] Christian | Multiple Time Scale Dynamics[END_REF] Sect. 5.4).

Center manifold reduction of the slow flow at the left fold point

The center manifold reduction technique (as presented for example in [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF][START_REF] Berlund | Noise-Induced Phenomena in Slow-Fast Dynamical Systems[END_REF]) is used in this section to reduce the slow flow [START_REF] Bergeot | Dynamic behavior analysis of a mechanical system with two unstable modes coupled to a single nonlinear energy sink[END_REF], in the vicinity of the left fold point (r LF , s LF , ∆ LF ), to the normal form of the dynamic saddle-node bifurcation which is then solve to deduce the nontrivial scaling law previously mentioned.

First, assuming = 0 in the right-hand sides of Eqs. (24b) and (24c), one considers the following simplification of the slow flow

r = f (r, s, ∆) (38a) s = g1 (r, s, ∆, 0) = g 1 (r, s, ∆), (38b)
∆ = g1 (r, s, ∆, 0) = g 2 (r, s, ∆). ( 38c
)
The goal in using the center manifold theorem is to reduce the dynamics with respect to the fast variables s and ∆, i.e. obtain an equivalent system which only depends on one variable. To achieve that, the Jacobian matrix J g of (38b-38c) evaluated at the fold point (r LF , s LF , ∆ LF ) is computed as

J g (r LF , s LF , ∆ LF ) =       -µ 2 1 9 1 -3µ 2 + 1 2- √ 1-3µ 2 α - 3 √ 1-3µ 2 -1 4 2- √ 1-3µ 2 α -µ 2       . ( 39
)
which has two real eigenvalues λ a = 0 and λ b = -µ.

Then to continue one must obtain polynomial nonlinearities. Indeed, the center manifold reduction requires isolating the linear terms. The functions g 1 and g 2 are therefore Taylor expanded around the fold point (r LF , s LF , ∆ LF ), up to order 2 for the variables s and ∆ and up to order 1 for the variable r. Remembering that by definition g 1 (r LF , s LF , ∆ LF ) = 0 and g 2 (r LF , s LF , ∆ LF ) = 0 that yields

g 1 (r, s, ∆) ≈ s -s LF ∂ s g 1 (r LF , s LF , ∆ LF ) + ∆ -∆ LF ∂ ∆ g 1 (r LF , s LF , ∆ LF ) + G1 (r, s, ∆) (40a) g 2 (r, s, ∆) ≈ s -s LF ∂ s g 2 (r LF , s LF , ∆ LF ) + ∆ -∆ LF ∂ ∆ g 2 (r LF , s LF , ∆ LF ) + G2 (r, s, ∆). (40b)
The polynomial functions G1 and G2 contain nonlinear terms with respect to the (s -s LF ) and (∆ -∆ LF ) (only quadratics terms are taken) and linear terms in (r -r LF ).

Denoting t = (s, ∆) T and t LF = (s LF , ∆ LF ) T and using Eq. ( 40), Eqs. ( 38b) and (38c) can be written in matrix form as

t = t -t LF J g (r LF , t LF ) + G r -r LF , t -t LF (41)
where G = ( G1 , G2 ) T . Now, to diagonalize the linear part of Eq. ( 41) the biorthogonal transformation is used again4 . For that, the right and left eigenvector matrices R = [r a r b ] and L = [l a l b ] are considered corresponding to the following right and left eigenvalue problems with respect to

J g (r LF , t LF ) J g (r LF , t LF )R = RΛ and J g (r LF , t LF ) T L = LΛ. ( 42
)
where Λ = diag(λ a , λ b ) = diag(0, -µ) denotes the diagonal matrix of the eigenvalues. Again the biorthogonal transformation consists of stating

t -t LF = Rq with q = L T t -t LF (43)
where q = (q a , q b ) T . Note that here, on a system undergoing a saddle-node bifurcation, the eigenvalues, the eigenvectors and the modal coordinates are real.

The variable q a is associated to λ a = 0 and the variable q b is associated to λ b = -µ. Moreover, the matrices R and L have been chosen such that L T R = I and L T J g (r LF , t LF )R = Λ. Eq. ( 41) is then expressed with respect to the variables q a and q b as

q = Λq + G (u, q) (44)
where G (u, q) = L T G r -r LF , Rq and u = r -r LF . Because of the latter change of variable the linear terms in q a and q b vanish in G.

The Jacobian matrix J g (r LF , t LF ) has a null eigenvalue and a negative one. Consequently, the center manifold theorem states that, in the neighborhood of the fold point (r LF , s LF , ∆ LF ), studying the following reduced system

q a = G 1 (q a , (q a ), u) (45)
is equivalent to studying Eq. ( 44). The function cannot be obtained explicitly and it must therefore be approximated. For sake of simplicity one chooses the tangent space approximation, i.e.

q b = (q a ) = 0.
To obtain a normal form of the dynamic saddle-node bifurcation final approximations are needed. First, because, we study the dynamics near the fold point (r LF , s LF , ∆ LF ) we assume that f (r, s, ∆) = f (r LF , s LF , ∆ LF ) = f LF . Then, only terms of order 2 in q a and of order 1 in u are kept in G 1 leading to

G 1 ≈ a 1 q 2 a + a 2 u (46)
with

a 1 = 3 1 -3µ 2 -1 µ 2 + 1 -3µ 2 + 1 18µ 2 and a 2 = 3α 2- √ 1-3µ 2 α 4 √ 2 3µ 2 + 1 -3µ 2 + 1 . ( 47
)
An approximated formulation of the slow flow (38) as the normal form of the dynamic saddlenode bifurcation is then finally obtained as

ˆ q a = h(q a , v) = q 2 a + v, ( 48a
) v = 1 (48b) with v = a 2 a 1 u, ˆ = f LF a 2 a 2 1 . ( 49
)
and where h(q a , v) satisfies the bifurcation conditions h(0, 0) = 0 and ∂ qa h(0, 0) = 0 and the conditions of the saddle-node bifurcation ∂ qaqa h(0, 0) = 0 and ∂ v f (0, 0) = 0. The bifurcation is said to be dynamic because unlike the static saddle-node bifurcation the bifurcation parameter v is a slowly varying parameter (because of (48b)).

The critical manifold M 0 of (48) consists of an attracting branch

M 0,a = {(x, y) ∈ R -× R -| x = - √ -y} and a repelling branch M 0,r = {(x, y) ∈ R + × R -| x = √ -y}.
This critical manifold is hyperbolic except in (0, 0). Finally, it may be noted that the center manifold reduction amounts to stating that locally the critical manifold is approximated by a parabolic function.

Analytical solution for the normal form of the dynamic saddle-node bifurcation

The method to solve analytically Eq. ( 48) is presented for exemple in [START_REF] Berlund | Noise-Induced Phenomena in Slow-Fast Dynamical Systems[END_REF]. It begins by dividing Eq. (48a) by Eq. (48b), that leads to

ˆ dq a dv = q 2 a + v. ( 50
)
The subsequent scaling q a = ˆ 1/3 z and v = -ˆ 2/3 s is then introduced to obtain the following equation

dz ds = -z 2 + s. ( 51
)
Setting z(s) = ϕ (s)/ϕ(s) in (51) yields the linear second order equation ϕ (s) = sϕ(s), whose solution can be expressed in terms of Airy functions as ϕ(s) = C 1 Ai(s) + C 2 Bi(s) where C 1 and C 2 are constants. Assuming that x(-∞) = -√ -y, only the contribution of Ai(s) is kept (this is not proved here). Therefore the solution with respect to the original variable q a , denoted q a (y), is then given by

q a (v) = ˆ 1/3 Ai -ˆ -2/3 v Ai -ˆ -2/3 v . ( 52
)
Note that from Eq. (48) to Eq. ( 52) no asymptotic approach has been used.

B. Bergeot

The previous analysis reveals the nontrivial -dependence of the slow flow near the fold point, involving the fractional exponents 1/3 and 2/3. Consequently, near the fold point, the trajectory no longer tracks the attracting branch of M 0,a of the critical manifold at a distance of order as it is the case when the trajectory evolves slowly far from the fold point (as predicted by the Fenichel theory in the normally hyperbolic branches of M 0 ).

Two particular values of r can be deduced from Eq. (55): the value at which s (r) = s LF , denoted r 0 , and the value for which s (r) tends toward infinity, denoted r ∞ . These values of r correspond to the first zeros of Ai and Ai respectively (see Eqs. ( 53) and ( 54)). Therefore, one has

s (r) = s LF for r = r 0 = r LF + K 0 2/3 (56)
where

K 0 = 1.01879 × a 1 a 2 f LF a 2 a 2 1 2/3 and s (r) → ∞ for r = r ∞ = r LF + K ∞ 2/3 (57)
where

K ∞ = 2.33810 × a 1 a 2 f LF a 2 a 2 1 2/3 .

New theoretical prediction of the mitigation limit

From the results obtained in the previous section, especially Eq. ( 57), it is possible to obtain a more accurate prediction of the mitigation limit than (36) which takes into account the value of the small parameter .

As seen previously (see Section 3.5), the mitigation limit is the value of ρ for which the fixed point s e 2 becomes smaller than s A (the abscissa of the arrival point of the trajectory on the critical manifold) or for which s e 2 no longer exists (see Eq. ( 34)). Within the zeroth-order approximation the arrival point is s A = s U defined by [START_REF] Berlund | Noise-Induced Phenomena in Slow-Fast Dynamical Systems[END_REF].

It is not easy to derive a new expression of s A from the scaling law (55). However, one can assume that when the trajectory reaches the right attracting par of M 0 , the limit value r ∞ has already been reached. Therefore, the solution of s e 2 = s A , denoted is equivalent to ρ U , is also the solution of

H s e 2 (ρ U ) = r ∞ , (58a) = r LF + K ∞ (ρ U ) 2/3 , (58b) = H(s LF ) + K ∞ (ρ U ) 2/3 . ( 58c
)
Except for very small values of for which ρ U ≈ ρ U 0 and the zeroth-order approximation is sufficient, the special value of µ solution of ρ U = ρ S is in general larger than 1/ √ 3, the upper limit of the study interval (see Eq. ( 30)). Therefore, we consider here that the mitigation limit is directly ρ U . From Eq. ( 58) a new and more accurate value of the optimal value of the NES damping coefficient can be derived. Indeed, denoting

H (s e 2 (ρ)) = f 1 (ρ, µ), H(s LF ) = f 2 (µ) and K ∞ (ρ) 2/3 = f 3 (ρ, µ) one has dρ µ = d µ f 2 + ∂ µ f 3 -∂ µ f 1 ∂ µ f 1 -∂ µ f 3 . ( 59
)
Therefore solving d µ ρ = 0 together with Eq. (58c) gives the coordinates of the maxima of ρ U (µ) and then optimal value of the NES damping coefficient, denoted µ opt .

Application to an aeroelastic aircraft wing system coupled to an NES

The methodology presented in previous sections is illustrated here on an aeroelastic aircraft wing system coupled to an NES (see Fig. 3) which has already been used in the past in several works about LCO mitigation by means of one or several NESs [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][START_REF] Bergeot | Steady-state regimes prediction of a multi-degree-of-freedom unstable dynamical system coupled to a set of nonlinear energy sinks[END_REF]. First, equations of motion the wing-NES system are presented in Section 5.1. Then, the comparison between theoretical results and numerical simulations is performed in Section 5.2.

Equations of the aircraft wing model coupled to one NES

The equations of motion of the model presented in Fig. 3 can be derived for example by means of the Lagrange equations. Assuming a small angle between the wing and the horizontal, the Lagrange equations yield the following equations of motion with respect to the heave z and the angle of attack (pitch) φ of the wing and to the NES displacement y (see Fig. 3)

mz + S ϕ φ + K z z + C z ż + qP a Γ ( φ + ż/U ) + K NL z z 3 + C y ż -d φ -ẏ + K NL y (z -d φ -y) 3 = 0 ( 66a
)
I ϕ φ + S ϕ z + K ϕ φ + C ϕ φ -qeP a Γ ( φ + ż/U ) + K NL ϕ φ3 + dC y ż -d φ -ẏ + dK NL y (z -d φ -y) 3 = 0 (66b) m y ÿ + C y d φ + ẏ -ż + K NL y (d φ + y -z) 3 = 0 (66c)
where again the time derivative is denoted " . ". b = L/2 is the semichord length. A is the aerodynamic center, B the elastic axis, G the center of gravity of the aircraft wing. e is the location aerodynamic center A measured from B (positive ahead of B). The quantities m and I ϕ are respectively the mass of the wing and its moment of inertia with respect to B. S ϕ = mBG is the mass unbalance in the wing. K z and K ϕ are the linear heave and pitch stiffnesses respectively whereas K NL z and K NL ϕ are the cubic heave and pitch stiffnesses. C z and C ϕ are the heave and pitch damping coefficients. U is the constant and uniform flow speed around the wing and q = 1 2 ρ ∞ U 2 is the dynamic pressure where ρ ∞ is the density of the flow. P a is the planform of the wing, Γ is the lift curve slope and d is the offset attachment of the NES to the wing. Finally, m y , C y and K NL y (m = 1, . . . , M ) are the mass, the damping coefficient and the cubic stiffness of the NES respectively.

For convenience, Eq. ( 66) is written in a dimensionless form as follows

x + s ϕ φ + Ω 2 x + ζ x x + γΘ Θ φ + x + ξx x3 + mh h = 0 (67a) r 2 ϕ φ + s ϕ x + r 2 ϕ φ + ζ ϕ φ -ηγΘ Θϕ + x + ξϕ φ3 + δ mh h = 0 (67b) mh h + ζh δ φ + h -x + ξh (δ φ + h -x) 3 = 0 (67c)
where the time t has been replaced by the dimensionless time t = ω ϕ t (with

ω ϕ = K ϕ /I ϕ
x + s ϕ ϕ + Ω 2 x + ζ x x + γΘ Θϕ + x + m h h = 0 (69a) r 2 ϕ ϕ + s ϕ x + r 2 ϕ ϕ + ζ ϕ ϕ -ηγΘ Θϕ + x + δm h h = 0 (69b) h + ζ h δϕ + h -x + ξ h (δϕ + h -x) 3 = 0 (69c)
with ζ h = ζh /m h and ξ h = ξh /m h . The physical bifurcation parameter under consideration is the reduced speed of the air flow Θ and the following parameters for the primary system are used ζ x = 0.01, ζ ϕ = 0.01, s ϕ = 0.2, r ϕ = 0.5, Ω = 0.5, η = 0.4, and γ = 0.2.

(

) 70
With the parameters (70), the aircraft wing model undergoes an Hopf bifurcation at Θ wo Hopf = 0.933 at which the real part of the first eigenvalue of Eq. ( 69) vanishes, i.e. Re(λ 1 ) = 0, and then becomes positive. It is recalled here that the pulsation ω and the generalized bifurcation parameter ρ are defined from the eigenvalue λ 1 of the unstable mode by Eq. ( 14).

The NES setup is defined by

m h = 1, δ = -0.8, ξ h = 5 (71)
and ζ h and which are specified below depending of the considered example. The several response regimes described in Section 3.4 are illustrated here by means of numerical simulations of the wing-NES system (67) (see Fig. 4) with, in addition to (70) and ( 71 4(d)) corresponding respectively to complete suppression, mitigation through a periodic mitigation through a SMR and finally no mitigation of the harmful LCOs of the primary structure. For comparison purposes, for each situation, a direct numerical integration of the aircraft wing system without NES is also shown.

Comparison between theoretical results and numerical simulations

The slow flow dynamics of Eq. (69), its critical manifold M 0 [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF], the zeroth-order analytical expressions of the mitigation limit (36) and of the optimal value of NES damping coefficient (37) are derived following the method described in Sections 3.1, 3.2 and 3.5 respectively. Note that through Eq. ( 19), the parameters µ and α are here defined as µ = ζ h /ω and α = ξ h /ω 2 . Simultaneously, the scaling law of the trajectory in the neighborhood of the left fold point of M 0 , i.e. Eq. ( 55), is computed. Afterwards, solving Eq. ( 58) and finding the root of Eq. ( 59), ρ U (µ) and µ opt are obtained. Finally, the asymptotic analytical expressions of the latter, given by Eqs. (64) and (65) respectively, are computed.

To illustrate the proposed analytical procedure, Fig. 5(a) shows the critical manifold M 0 (27), superimposed to the direct numerical integration of the slow flow (38) (see also Fig. 5(b) for the times series of s(τ ) and r(τ ) focused on the first relaxation cycle) and its scaling law given by Eq. ( 55). The values of r LF = H(s LF ) (see Eq. ( 29)) and r ∞ = r LF + K ∞ 2/3 (see Eq. ( 57)) are depicted by horizontal lines. The numerical integration of the slow flow (see red dashed line in Fig. 5(a)) shows a relaxation oscillations scenario in (s, r)-plane. Indeed, the trajectory of undergoes a succession The figure reveals a good agreement between the numerical simulation of the slow flow and its analytical scaling law in the neighborhood of the left fold point. One can see also that choosing r ∞ as the ordinate of point of intersection between the trajectory of the slow flow and the right attracting branch of M 0 is a more accurate approximation than choosing r LF as this is the case within the zeroth-order approximation.

In Fig. 6, the comparison in now performed in terms of mitigation limit. The theoretical predictions, the one obtained within the zeroth-order approximation (36) and those derived from the scaling law (solving numerically Eq. ( 58) and the one given by Eq. ( 64)) are compared with mitigation limits measured on direct numerical integration of the approximated slow flow (38) and on the full order system (67). In these last two cases, direct numerical integrations are performed for increasing values of ρ (by increasing Θ) and the mitigation limit corresponds to the last value of ρ allowing a harmless situation. In general one observes at this value the transition from a SMR to a "no mitigation regime" as shown for example in Fig. 4. All the mitigations limits are plotted with respect to the NES damping coefficient µ (with µ ∈ [0, 1/ √ 3] 5 ) and for four values of the perturbation parameter . One takes first a very little value ( = 0.001) to be sure to respect the assumptions of the asymptotic analysis (i.e. 0 < 1) and thus validate the method under these conditions. Then, three larger values are taken ( = 0.005, 0.02 and 0.1) to evaluate the robustness of the method when one deviates from the assumption of a small . To relate these results to the physical bifurcation parameter Θ, the latter is plotted with respect to the generalized bifurcation parameter ρ in Fig. 7 for the same values of as in Fig. 6. First, in Fig. 6, one can see that the zeroth-order approximation gives an inaccurate prediction, even for the smallest value of = 0.001. Moreover, one observes that, on the contrary, solving Eq. (58) provides a very good prediction of the mitigation limits measured on numerical simulations. For the two following values of (0.005 and 0.02), these theoretical prediction is closer to the mitigation limits obtained from numerical simulations of the full order system (67) than to those measured on numerical integrations of the approximated slow flow (38). This is surprising because Eq. ( 58) is derived from Eq. (38) but, in general, in these kind of unstable systems coupled to NES, the mitigation limits measured on the full order system are slightly smaller than those measured on the slow flow. Here, when increases the prediction depreciates becoming smaller than mitigation limits obtained from numerical simulations of Eq. ( 38) and getting closer to those measured on full order system (67). For the largest value = 0.1, the theoretical mitigation limit obtained solving Eq. (58) (in green on the figure) depreciates again but it is still satisfactory in a qualitative point of view. Moreover one can see that the theoretical prediction (in green) no longer has a local maximum and the optimal value of µ is equal to 1/ √ 3, its limit value for the existence of the SMR regimes (see Eq. ( 30)).

Finally, the prediction provided by Eq. ( 64) deteriorates enormously for large values of . This may be due to the non-convergent nature of the series (61), as already mentioned above (see comments after Eq. ( 64)).

Finally, in Fig. 8, the theoretical optimal value of the NES damping coefficient obtained from the zeroth-order approximation µ opt 0 (see Eq. ( 37)) is compared to those derived from the scaling law (55): the one obtained solving numerically d µ ρ = 0 (with d µ ρ given by Eq. ( 59)) together with Eq. (58c) and the asymptotic expression given by Eq. (65). All this is plotted with respect to the perturbation parameter . Assuming a good agreement between the theoretical prediction of the mitigation limit obtained solving Eq. (58) and the actual value (as shown by Fig. 6), only theoretical values are presented in Fig. 8. The figure highlights the relevance of the proposed results especially for large values of the parameter . We can also see that the asymptotic expression is here able to describe qualitatively the evolution of µ opt with respect to the perturbation parameter for 0.04. For 0.04, the prediction obtained from Eqs. (58c) and (59) tends to 1/ √ 3 whereas the prediction given by Eq. (65) continues to grow.

Conclusion

In this work a mechanical system with one unstable mode coupled to a Nonlinear Energy Sink (NES) has been studied. The equations of motion of the system have been first simplified using a reducedorder model for the primary structure by keeping only its unstable modal coordinates. Introducing a small perturbation parameter related to the mass ratio between the NES and the primary structure, the slow flow has been then derived by means of the complexification-averaging method. Because of the presence of the small perturbation parameter, the slow flow is governed by two different time scales. More precisely, within the framework of the geometric singular perturbation, in its real form it appears as a (2, 1)-fast-slow system. The critical manifold of this slow flow has been computed. Within the zeroth-order approximation, i.e. the analysis within the limit case in which the perturbation parameter is equal to zero, the solution of the slow flow is approximated by a discontinuous phase trajectory in which the fast and slow epochs are described independently of each other. In particular, the slow epoch are supposed to occur on attracting parts of the critical manifold. Moreover, in this context a theoretical prediction of the mitigation has been obtained. The mitigation limit being defined as the value of the generalized bifurcation parameter, which is obtained rescaling the eigenvalue real part of the primary structure corresponding to the unstable mode, which separates harmless situations (in which the NES acts) from the harmful situation (in which the NES does not act). The results obtained within the zeroth-order approximation (both phase trajectory and mitigation limit) depreciate for the largest values of the perturbation parameter (the latter remaining small in the context of this work).

The novelty of the present work is to use the center manifold theorem to reduce the slow flow, in the neighborhood of a fold point of its critical manifold, to a normal form of the dynamic saddlenode bifurcation. This allowed us to obtain a scaling law for the slow flow in this neighborhood, i.e. a law which describes the dependance, with respect to the small perturbation parameter, of the distance between the critical manifold and the actual trajectory of the slow flow. The law reveals given by Eq. (37) (blue line) is compared to those derived from the scaling law (55): the one obtained solving numerically d µ ρ = 0 (with d µ ρ given by Eq. ( 59)) together with Eq. (58c) (green line) and the asymptotic expression given by Eq. (65) (orange line). As previously, parameters (70) and (71) are used.

that the slow flow scales in a nontrivial way with respect to the perturbation parameter, involving the fractional exponents 1/3 and 2/3. Finally, a theoretical prediction of the mitigation limit has been deduced from the scaling law. With certain assumptions, it was shown that this prediction is a correction with respect to the prediction obtained within the zeroth-order approximations, involving also a fractional order dependence with respect to the perturbation parameter.

The last section of the paper has been devoted to illustrate and validate the proposed methodology on a an aeroelastic aircraft wing model coupled to one NES. The accuracy of the present theoretical results compared to those obtained on numerical simulations of the full order system together with the general nature of the method suggest that the latter may constitute a tool for the practical design of NES.

Figure 1 .

 1 Figure 1. Typical example of the critical manifold in the (s, r)-plane given by Eq. (28a) for µ = 0.25 and α = 5.

  ) and time derivative d t is denoted " ". The dimensionless displacements are defined by x = z/b, h = y/b. Moreover, δ = d/b, s ϕ = S ϕ /(mb), Ω = ω z /ω ϕ (with ω z = K z /m) and η = e/b. r ϕ = I ϕ /(mb 2 ) is the radius of gyration of the cross section of the wing. The dimensionless nonlinear stiffnesses and damping coefficients are ξx= K NL z b 2 /(mω 2 ϕ ), ζ x = C z /(mω ϕ ), ξϕ = K NL ϕ /(mb 2 ω 2 ϕ ), ζ ϕ = C ϕ /(mb 2 ω ϕ ), ξh = K NL y b 2 /(mω2 ϕ ), ζh = C y /(mω ϕ ). The mass and density ratios are defined B. Bergeot h = h √ to obtain the final form of the equations of motion

  ), = 0.01, ζ h = 0.3, ξ x = 4 and ξ ϕ = 8. Four values of the reduced speed are used, namely: Θ = 0.94 (see Fig. 4(a)), Θ = 0.945 (see Fig. 4(b)), Θ = 0.97 (see Fig. 4(c)) and Θ = 0.98 (see Fig.

Figure 5 .

 5 Figure 5. Illustration of the proposed analytical procedure. (a) The critical manifold M 0 (27) (black line) superimposed to the direct numerical integration of the slow flow (38) (red line) and its scaling law near the left fold point given by Eq. (55) (blue line, the dashed parts are the horizontal asymptotes of s (r) due to the zeros of the Airy function). The values of r LF = H(s LF ) (see Eq. (29)) and r ∞ = r LF + K ∞ 2/3 (see Eq. (57)) are depicted by magenta and green horizontal lines respectively. (b) Times series of s(τ ) (blue line) and r(τ ) (green line) obtained from the numerical integration of the slow flow (38). The times series are focused on the first relaxation cycle. Parameters (70) and (71) are used with ζ h = 0.25, = 0.005 and Θ = 0.95.

Figure 6 .

 6 Figure 6. Theoretical mitigation limits: the one obtained within the zeroth-order approximation (36) (blue line) and those derived from the scaling law (solving numerically Eq. (58) (green line) and the one given by Eq. (64) (orange line)) compared with mitigation limits measured on direct numerical integration of the approximated slow flow (38) (gray dots) and on the full order system (67) (black dots). Parameters (70) and (71) are used and ((a) = 0.001, (b), = 0.005, = 0.02 and (d) = 0.1.

Figure 7 .

 7 Figure 7. The physical bifurcation parameter Θ with respect to the generalized bifurcation parameter ρ. As in Fig. 6, parameters (70) are used and (a) = 0.001, (b), = 0.005, = 0.02 and (d) = 0.1.

Figure 8 .

 8 Figure8. The theoretical optimal value of the NES damping coefficient with respect to the small perturbation parameter . The one obtained from the zeroth-order approximation µ opt 0

Sometimes the terms "slow" and "super-slow"[START_REF] Gendelman | Bifurcations of self-excitation 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] are used to leave the term "fast" for the oscillations on which the average was carried out and therefore to be coherent with term slow flow. In the present work the terms "slow" and "fast" are preferred to be in agreement with the vocabulary usually used in the literature on dynamical systems.

The difference between critical and invariant manifolds will be calrified in the manuscript in Section 4.

These terms do not affect qualitatively the dynamic behavior of the slow flow.

To avoid having to introduce new notations, same symbols as in Section 2.2 are used in this section to describe the biorthogonal transformation except for the subscripts which are now letters (a, b,...) instead of numbers (1,2,...).

Remembering that in this work only situations in which Eq. (30) holds are considered.

Acknowledgments

We thank the anonymous reviewers for their careful reading of our manuscript and their many relevant comments and suggestions.

2.

Illustration of the normal form of the dynamic saddle-node bifurcation. Result of the numerical integration of Eq. ( 48) with initial condition (q a (0) = -1, v(0) = -0.5) (red dashed line) compared to the analytical scaling law q a (y) given by (52) (blue line, the dashed parts are the horizontal asymptotes of q a (y) corresponding to the zeros Airy function) for = 0.01. The first zero and the first singularity of q a (v) (orange and green lines respectively) and the corresponding critical manifold M 0 (attracting part M 0,a in black and repelling part M 0,r in gray) are also shown.

The first zero of q a (v) is given by the first zero of the Airy function derivative Ai and the first singularity by the first zero of Ai which are tabulated values (see e.g. [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th Edition, Applied mathematics series[END_REF]). Therefore, one has

and

An illustration of the normal form of the dynamic saddle-node bifurcation is shown in Fig. 2. The direct numerical integration of Eq. (48) with initial condition (q a (0) = -1, v(0) -0.5) is compared to the analytical solution q a (y) given by (52) for = 0.01. The first zero, the first singularity of q a (v) and the corresponding critical manifold M 0 (which consists of an attracting part M 0,a and a repelling part M 0,r ) are also shown.

From Eqs. ( 43) and (49), an approximated analytical solution of the slow flow (38) in the neighborhood of the fold point (r LF , s LF , ∆ LF ) is obtained as

where a 1 , a 2 and ˆ are given by ( 47) and (49) respectively and r a,1 is the first coordinate of the vector r a .

Introducing a new perturbation parameter ν = 2/3 , Eq. ( 58) takes the following form

Because 0 < 1 one has also 0 < ν 1, the latter equation is a perturbation of H (s e 2 ) = H(s LF ) which is equivalent to H (s e

2 ) = H(s U ) and consequently to s e 2 = s U whose solution is ρ U 0 given by Eq. ( 35).

Therefore, to find an explicit formulae for the solution of (58), a regular perturbation approach is used by expressing ρ U as a power series of ν

Then, substituting (61) into (60) and equating the coefficients of identical power of ν on both sides of the equation lead to the following set of equations (up to the first order)

As expected, from (62a) one has ρ U ,0 = ρ U 0 (see Eq. ( 35)) and therefore Eq. (62b) yields

Finally, to a first approximation, one has the following asymptotic expression for ρ U

The asymptotic series (61) converges if the functions of ρ in Eq. (60) (i.e. s e 2 (ρ) and K ∞ (ρ)) can be expanded as a convergent Taylor series at ρ = ρ U 0 . Unfortunately s e 2 (ρ) has an infinite derivative at ρ = ρ S (see Eq. ( 34)). Because of that, at µ = µ * where ρ U 0 = ρ S the term ρ U ,1 in the series Eq. ( 61) is equal to zero but the following terms grow to infinity making the series divergent. However, except at µ = µ * and for sufficiently small values of , the expression (64) shows again a nontrivial -dependence stating that ρ U is a O( 2/3 )-perturbation of ρ U 0 . Equating to zero the derivative, with respect to µ, of the right-hand side of Eq. ( 64) and performing again a perturbative technique with respect to ν = 2/3 , an approximative expression of µ opt is obtained as

where µ opt 0 is given by Eq. (37). Obviously, the nontrivial -dependence is found again as, if is sufficiently small, µ opt is a O( 2/3 )-perturbation of µ opt 0 . Of course the limitations previously mentioned for Eq. ( 64) are also valuable for Eq. (65). by mh = m y /(mω 2 ϕ ) and γ = bP a ρ ∞ Γ/(2m) respectively. Finally, the bifurcation parameter under consideration is the reduced speed of the flow Θ = U/(bω ϕ ). Equation (67) has the general form given by Eq. ( 2) with x = (x, φ) T and .98 corresponding respectively to complete suppression, mitigation through a periodic regime, mitigation through a SMR and finally no mitigation of the harmful limit cycle oscillations of the primary structure. For comparison purposes, for each situation, a direct numerical integration of the wing system without NES is also shown (in green).

Zoom into the NES installation

of almost horizontal fast parts and slow parts evolving in the O( )-vicinity of the critical manifold M 0 (here, at the scale of the figure and apart from the fold points, the trajectory and M 0 seem superimposed). The first horizontal fast part is from the initial condition (s(0) = 0.15, r(0) = 0.04) to the left attracting part of M 0 and the other are fast jumps between the left and right attracting