with $A_0 = A_0^R + jA_0^I$, $A_1 = A_1^R + jA_1^I$ and $A_{1m} = A_{1m}^R + jA_{1m}^I$ where $A_{1m}$ is the m-component of the vector $\mathbf{A}_1$.

Hence, finding the fixed points of Eq. (A.35) goes back to find the roots of

$$ (A_0^R + A_1^R H_1(r_1)^2) H_1(r_1) + \sum_m \sum_m^M r_m^2 (A_{1m}^R R_m(r_m) + A_{1m}^I I_m(r_m)) = 0 \quad (A.66a) $$

$$ H_m(r_m) - H_1(r_1) = 0, \quad (m = 2, \dots, M), \quad (A.66b) $$

which can be easily solved with a regular computer (at least until $N \approx 5$ or 6).

To check the stability of a fixed point (denoted $\mathbf{r}^* = (r_1^*, \dots, r_N^*)$), Eq. (B.40) is written as

$$ \frac{\partial r_1}{\partial \tau} = f(r_1, \dots, r_N), \quad (A.67) $$

where

$$ f(r_1, \dots, r_N) = \frac{2 (A_0^R + A_1^R H_1(r_1)^2) H_1(r_1) + 2 \sum_m^M r_m^2 (A_{1m}^R R_m(r_m) + A_{1m}^I I_m(r_m))}{\frac{dH_1}{dx}(r_1)}, \quad (A.68) $$

and the stability is deduced from the sign of $\frac{df}{dr_1}|_{\mathbf{r}=\mathbf{r}^*}$ where, using Eq. (A.66b), $\frac{df}{dr_1}|_{\mathbf{r}}$ takes the form

$$ \left. \frac{df}{dr_1} \right|_{\mathbf{r}} = \frac{\partial f}{\partial r_1}(\mathbf{r}) + \sum_{m=2}^{M} \frac{\partial f}{\partial r_m}(\mathbf{r}) \frac{dr_m}{dr_1} \quad (A.69a) $$

$$ = \frac{\partial f}{\partial r_1}(\mathbf{r}) + \sum_{n=2}^{M} \frac{\partial f}{\partial r_m}(\mathbf{r}) \frac{\frac{dH_1}{dx}(r_1)}{\frac{dH_m}{dx}(r_m)}. \quad (A.69b) $$

Moreover, to be stable a fixed point must be on a stable branch of the critical manifold S.

Finally, a fixed point $\mathbf{r}^* = (r_1^*, \dots, r_M^*)$ of Eq. (B.40) is **stable** if the two following conditions are satisfied

$$ 1. \quad \left. \frac{df}{dr_1} \right|_{\mathbf{r}=\mathbf{r}^*} < 0, $$

$$ 2. \quad \forall m, \ 1 \le m \le M, \ r_m^* \in [0, r_m^{max}] \cup [r_m^{min} + \infty). $$

A stable (respectively unstable) fixed point is denoted $\mathbf{r}_s^* = (r_{s,1}^*, \dots, r_{s,M}^*)$ (respectively $\mathbf{r}_u^* = (r_{u,1}^*, \dots, r_{u,M}^*)$).

Finally, each fixed point (stable or unstable) as a point of S is a part of one of branches $S_k$.

### A.3.5 Prediction of the steady-state regimes

As discussed for example in [a7, a14, 75, 78, 123] when a NES is coupled to an unstable system, four types of steady-state regimes can be generated corresponding to a complete suppression of the instability, mitigation through Periodic Response (PR), mitigation through Strongly Modulated Response (SMR) or no mitigation. Complete suppression means that the trivial fixed point, common to both the non-averaged system and the slow-flow, is stable. Mitigation through PR means that a stable (nontrivial) fixed point of the slow-flow is reached leading, for non-averaged system, to a LCO with an amplitude smaller than the amplitude of the corresponding LCO undergone by the primary system without absorbers. Mitigation through SMR means that the slow-flow undergoes relaxation oscillations corresponding to a phase-amplitude modulated regime for the non-averaged system, the so-called "Strongly modulated response". Finally no mitigation means that 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 with respect to the variable s. This stable fixed point denoted $\mathbf{r}_s^*$ corresponds to the stable fixed point at which the amplitude with respect to the variable s is maximized. This stable fixed point is positioned in the branch $S_{2M}$.

As proposed in [a7], these regimes can also be classified into two groups named harmless situation and harmful situation reflecting the situations in which the NESs act or not. Harmless situation includes complete suppression, mitigation through PR and mitigation through SMR regimes. Harmful situation corresponds to no mitigation regime. The main objective of NESs is to put the system in harmless situations.

The nature of the steady-state regime depends on some features of the slow-flow. Firstly, the initial conditions from which one can know where the fast dynamics leads the system trajectory on the critical
manifold. Secondly, the slow subsystem (A.39), in particular its fixed points (position and stability), which guides the dynamics on $S$ at the slow time scale. In accordance with real word situations, the following prediction is limited to a set of initial conditions [$s(0), r_1(0), \dots, r_N(0)$] as a small perturbation of the trivial solution.

To predict the steady-state regimes, we have the CM $S$, the stable domain $D$ of $S$ in terms of the Cartesian products $I_k$ and the associated branches $S_k$ in the $(r_1, \dots, r_M)$-space with the corresponding bounds $\mathbf{b}_k^1$ and $\mathbf{b}_k^2$ (which are fold points) and the associated arrival points $\mathbf{a}_k^1$ and $\mathbf{a}_k^2$. We have also the stable fixed points $\mathbf{r}_s^*$ and unstable fixed points $\mathbf{r}_u^*$ which are positioned in the corresponding branches $S_k$. Note that the bound point $\mathbf{b}_1^1$ is the trivial equilibrium of the system.

Recalling the two subsystems (A.38) and (A.39), the trajectory of the system is a succession of slow (given by (A.38)) and fast (given by (A.39)) parts. The slow parts lead in $S$ whereas the fast parts are jump occurring out of $S$ between two points of $S$.

From a given initial condition, after a first fast transient response (a jump), the trajectory reaches the critical manifold $S$ on the branch $b_1 = S_1$ at the point $\mathbf{p}_1$.

We assume that we are at the end of the step $l$ of the procedure i.e on the branch $b_l = S_{k_l}$ at the (arrival) point $\mathbf{p}_l$.

The step $l+1$ (from $l$ to $l+1$) reads as follows.

The trajectory follows $b_l = S_{k_l}$ slowly starting from $\mathbf{p}_l$ and the following cases are possibilities :

**Case 1.** There are no fixed points in $S_{k_l}$, then, depending on the sign of the function $f(r_1, \dots, r_N)|_{\mathbf{r}=\mathbf{p}_l}$ (see Eq. (A.67)), $\mathbf{b}_{k_l}^1$ or $\mathbf{b}_{k_l}^2$ are reached, the trajectory leaves $S_{k_l}$ jumping to the arrival point $\mathbf{p}_{l+1} = \mathbf{a}_{k_l}^1$ (or $\mathbf{a}_{k_l}^2$) on an other branch $S_{\bar{k}}$ given $k_{l+1} = \bar{k}$ and $b_{l+1} = S_{\bar{k}}$.

**Case 2.** The arrival point $\mathbf{p}_l$ is between an unstable fixed point and $\mathbf{b}_{k_l}^1$ (resp. $\mathbf{b}_{k_l}^2$), then $\mathbf{b}_{k_l}^1$ (resp. $\mathbf{b}_{k_l}^2$) is reached, the trajectory leaves $S_{k_l}$ jumping to the arrival point $\mathbf{p}_{l+1} = \mathbf{a}_{k_l}^1$ (resp. $\mathbf{a}_{k_l}^2$) on an other branch $S_{\bar{k}}$ given $k_{l+1} = \bar{k}$ and $b_{l+1} = S_{\bar{k}}$.

**Case 3.** The arrival point $\mathbf{p}_l$ is between an stable fixed point and $\mathbf{b}_{k_l}^1$ (resp. $\mathbf{b}_{k_l}^2$), then the stable fixed point is reached and the procedure stops.

**Case 4.** The arrival point $\mathbf{p}_l$ is between an stable fixed point and unstable fixed point, then the stable fixed point is reached and the procedure stops.

At the end of step $l+1$, if the trajectory reaches a new branch $b_{l+1}$ at the point $\mathbf{p}_{l+1}$ a new step can be performed.

The procedure stops when one of the following conditions are met resulting in harmless or harmful situations :

* The system is in a **harmless situation** if

- During a step, the trajectory stops on a stable fixed point which differs from $\mathbf{r}_s^{*m} \in S_{2^M}$ resulting to a Mitigation through Periodic Response (Cases 3 or 4).

- At each step, the trajectory reaches a new branch. Due to the finite number of branch $S_k$, the trajectory will reach the same branch in the same condition resulting to a Mitigation through Strongly Modulated Response (SMR) (Cases 1 or 2).

* The system is in a **harmful situation** if the trajectory stops on the stable fixed point $\mathbf{r}_s^{*m} \in S_{2^M}$. Indeed, in this case, the reached fixed point has a large amplitude close to that of the case without NES.

Using the notion of harmless and harmful situation, the concept of **mitigation limit** of the system can be introduced. For a given parameter of the system named "bifurcation parameter", the mitigation limit is defined as the value of a chosen bifurcation parameter which separates harmful situation from harmless situation.

From a given set of parameters, the procedure described hereinbefore gives a theoretical prediction of the resulting steady-state regime and therefore allows to know if the system is a harmless or harmful situation. Consequently, the mitigation limit can be predicted theoretically as the first value of the chosen bifurcation parameter for which harmful situation holds.

## A.4 APPLICATIONS

Two self-excited mechanical systems are used in this section to check the relevance of the analytical procedure proposed in Sect. A.3.5. The first one is the well-known Hultén's model [102] which reproduces the typical dynamic behavior of friction systems. The second one is a 2-DOF nonlinear airfoil model. It has been already used by Lee et al in [123].
<figure><img src="image_66.png" /><figcaption>FIGURE A.3. (a) Hultèn's model with two NES. (b) Zoom on the NES1. (c) Zoom on the NES2.</figcaption></figure>

### A.4.1 Mode coupling instability mitigation in a friction system using two NESs

The system is shown in Fig. F.1. It is composed of a 2-DOFs Hultèn's model (the primary system) coupled to two NESs in ungrounded configuration. This model has already been considered in [a14] where it was shown that the use of NES appears to be an interesting way to control mode-coupling instability in braking systems. The objective here is to verify that the method presented here allows us to distinguish between strongly modulated responses and no mitigation responses, the method in [a14] being not able to make this distinction in all the cases.

Following Sect. D.2.1, the equations of motion of the system take the form (see B.6)

$$x_1'' + \zeta_{x_1} x_1' + x_1 - \mu \Omega^2 x_2 + \epsilon (\xi_{x_1} x_1^3 - \mu \xi_{x_2} x_2^3) + \epsilon m_{h_1} h_1'' = 0 \quad (A.70a)$$

$$x_2'' + \zeta_{x_2} x_2' + \Omega^2 x_2 + \mu x_1 + \epsilon (\mu \xi_{x_1} x_1^3 + \xi_{x_2} x_2^3) + \epsilon m_{h_2} h_2'' = 0 \quad (A.70b)$$

$$h_m'' - \gamma_{h_m} (x_m' - h_m') - \alpha_{h_m} (x_m - h_m)^3 = 0 \quad (m = 1, 2) \quad (A.70c)$$

where Eqs. (A.70) correspond to Eqs. (D.3). The corresponding Eqs. (D.5) and Eqs. (D.10) hereafter referred indistinctly as the Reference Model can be easily deduced as well as the slow-flow model corresponding to Eqs. (A.31).

The bifurcation parameter under consideration throughout this section is the friction coefficient $\mu$. The following set of numerical values are used, for the other parameters: $\xi_{x_1} = 5$, $\xi_{x_2} = 0$, $\zeta_{x_1} = 0.01$, $\zeta_{x_2} = 0.08$ and $\Omega = 0.9$ for the primary system, and $\gamma_{h_1} = 0.75$, $\gamma_{h_2} = 0.2$, $\alpha_{h_1} = 4$ and $\alpha_{h_2} = 9$ for the NESs. Moreover, $m_{h_1} = m_{h_2} = 2$ and several values for $\epsilon$ will be considered.

This set of numerical values leads to a Hopf bifurcation for the Hultèn's model without NESs at $\mu_{\text{Hopf}}^{\text{wo}} = 0.073$. If $\mu < \mu_{\text{Hopf}}^{\text{wo}}$ the system is stable and if $\mu > \mu_{\text{Hopf}}^{\text{wo}}$ the system is unstable. In the latter case, the model variables increase and reach a LCO due to the nonlinear terms. The LCOs are what we call "harmful situations", i.e. the large amplitude vibrations we want to mitigate using the two NESs. Here we focus on the possible mitigation regimes which can appear when the primary system is unstable.

First of all, the procedure described in Sect. A.3.5 for the prediction of the trajectory of the slow-flow in the $(r_1, r_2)$-space is illustrated in Fig. A.4 for $\epsilon = 2 \cdot 10^{-3}$ and four values of the bifurcation parameter $\mu$ (friction coefficient). The CM $S$, computed from (A.47) and (A.52), is depicted with their stable parts (in gray), i.e. the $S_k$ (with $k = 1, \dots, 2^M = 4$), and unstable parts (in black). The stable and unstable fixed points, computed from (B.43) and (A.69), are represented by magenta and blue points respectively. We can see in Figs. A.4(b) to A.4(d) a stable fixed point with a large value on the unbounded straight line (i.e. $S_{2^M} = S_4$). The latter corresponds to a large amplitude LCO mentioned above, if it is reached, we are in a harmful situations.

Following the procedure (as described in Sect. A.3.5) the theoretical trajectory of the slow-flow is determined and plotted with a dashed green line. This theoretical trajectory is compared to the trajectory obtained from the direct numerical integration of the slow-flow (plotted in red). The Fig. A.4(a) shows an harmless situation through a periodic regime, i.e. a stable fixed point of the slow-flow (in magenta) which is not in $S_4$.
<figure><img src="image_67.png" /><figcaption>FIGURE A.4. Critical manifold in the $(r_1, r_2)$-space defined through Eq. (E.32) together with theoretical and numerical trajectories of the slow-flow (A.31). Selected parameter values with $\epsilon = 2 \cdot 10^{-3}$.</figcaption></figure>

(here it is in $S_1$) is reached. In this case, the theoretical and numerical trajectories of the slow are almost superimposed. In Figs. A.4(b) and A.4(c) again harmless situations are depicted but through SMRs. Indeed, two different scenarios of relaxation oscillations of the slow-flow in $r_1$-direction are observed. Finally, because a stable fixed point is reached in $S_4$, the Fig. A.4(d) corresponds to an harmful situation. One can see that for the last three situations (i.e Figs. A.4(b), A.4(c) and A.4(d)), during the slow epochs on $S$ the theoretical and numerical trajectories of the slow-flow are again almost superimposed. However, during fast epochs (i.e. the jumps between each slow epoch) a difference is observed. In particular, unlike the theoretical trajectory, the numerical trajectory is not at each jump parallel to only one direction of the $(r_1, r_2)$-space. This difference is due to the 0-order approximation of the proposed asymptotic approach. Indeed, figures similar to Fig. A.4 using $\epsilon = 10^{-4}$ would show a perfect agreement between theoretical and numerical trajectories.

In addition, Fig. A.5 shows the direct numerical integration of the Reference Model (solid blue line) and of the slow-flow (red line) for situations corresponding to Figs. A.4(c) and A.4(d). In Fig. A.5(a), we can see that only $w_1(t)$ is a strongly amplitude modulated signal (i.e. SMR) which corresponds effectively to relaxation oscillations of the slow-flow in the $r_1$-direction. A large amplitude periodic regime (which is a harmful situation) corresponding to a large stable fixed point of the slow-flow is observed for the reference model signals in Fig. A.5(b).

Now, the relevance of the theoretical estimation of the trajectory with respect to the parameter $\epsilon$ is investigated. The comparison is first performed in term of the mitigation limit. For this purpose, the following quantities are defined:

* $\mu_{ml}^{th}$: the theoretical value of the mitigation limit,

* $\mu_{ml,SF}^{num}$: the numerical estimation of the mitigation limit "measured" on the graphs of the maximum steady-state amplitude obtained from numerical simulations of the slow-flow,

* $\mu_{ml,RM}^{num}$: the numerical estimation of the mitigation limit "measured" on the graphs of the maximum steady-state amplitude obtained from numerical simulations of the Reference Model.

Then the comparison is completed checking the differences in term of NES efficiency rate, denoted $\Lambda$ (in %) and defined as the relative difference between the mitigation limit and the Hopf bifurcation point of the airfoil model without NES

$$ \Lambda_{ml}^{th} = 100 \times \frac{\mu_{ml}^{th} - \mu_{Hopf}^{wo}}{\mu_{Hopf}^{wo}}, \qquad (A.71a) $$
<figure><img src="image_68.png" /><figcaption>FIGURE A.5. Direct numerical integration of the Reference Model (solid blue line) and of the slow-flow (red dashed line). Parameters: the same as those used for (a) Fig. A.4(c) and (b) Fig. A.4(d).</figcaption></figure>

$$ \Lambda_{\text{ml,SF}}^{\text{num}} = 100 \times \frac{\mu_{\text{ml,SF}}^{\text{num}} - \mu_{\text{Hopf}}^{\text{wo}}}{\mu_{\text{Hopf}}^{\text{wo}}} \qquad (\text{A.71b}) $$

$$ \Lambda_{\text{ml,RM}}^{\text{num}} = 100 \times \frac{\mu_{\text{ml,SRM}}^{\text{num}} - \mu_{\text{Hopf}}^{\text{wo}}}{\mu_{\text{Hopf}}^{\text{wo}}}. \qquad (\text{A.71c}) $$

This quantity gives indeed a direct estimation of the efficiency of the NESs: the larger it is, the more efficient the NESs network is.

The results of the comparison are reported Tab. A.1. One can observe that the mitigation limits obtained from numerical simulations of the slow-flow and of the reference model are very close. Moreover, the comparison between theoretical and numerical values shows a good agreement. Indeed, a difference of the order of magnitude of $\epsilon$ is observed. However, in terms of NES efficiency rate, we can see that in efficient situations (i.e. $\epsilon = 2 \cdot 10^{-2}$) the theoretical prediction can lead to a relatively large overestimation (a relative error of 16.28 %). As previously, these overestimations are due to the 0-order approximation of the asymptotic approach presented in Sect. A.3. This could be reduced by considering larger orders of approximation.

### A.4.2 Flutter mitigation in an airfoil system using M NESs

The system is shown in Fig. D.3. It is composed by a 2-DOF nonlinear airfoil model coupled to M NESs. Following Sect. D.2.1, the equations of motion of the system take the form (see A.7)

$$ x'' + s_{\varphi}\varphi'' + \Omega^2 x + \epsilon\xi_x x^3 + \zeta_x x' + \mu\Lambda\Theta(\Theta\varphi + x') + \epsilon \sum_{m=1}^{M} m_{h_m} h_m'' = 0 \qquad (\text{A.72a}) $$

$$ r_{\varphi}^2 \varphi'' + s_{\varphi} x'' + r_{\varphi}^2 \varphi + \epsilon \xi_{\varphi} \varphi^3 + \zeta_{\varphi} \varphi' - \eta \mu \Lambda \Theta (\Theta \varphi + x') + \epsilon \sum_{m=1}^{M} \delta_m m_{h_m} h_m'' = 0 \qquad (\text{A.72b}) $$

$$ h_m'' + \gamma_{h_m} (\delta_m \varphi' + h_m' - x') + \alpha_{h_m} (\delta_m \varphi + h_m - x)^3 = 0 \quad (m = 1, \dots, M) \qquad (\text{A.72c}) $$
<table><caption>**TABLEAU A.1.** Comparison between numerical and theoretical estimations of the mitigation limit and the NES efficiency rate for several values of the little parameter $\epsilon$.</caption><thead><tr><th></th><th></th><th colspan="3">ε</th></tr><tr><th></th><th></th><th>2 · 10<sup>-4</sup></th><th>2 · 10<sup>-3</sup></th><th>2 · 10<sup>-2</sup></th></tr></thead><tbody><tr><td rowspan="3">Mitigation limit</td><td>μ<sub>ml</sub><sup>th</sup></td><td>0.07553</td><td>0.0929</td><td>0.173</td></tr><tr><td>μ<sub>ml,SF</sub><sup>num</sup></td><td>0.07546</td><td>0.0918</td><td>0.162</td></tr><tr><td>μ<sub>ml,RM</sub><sup>num</sup></td><td>0.07543</td><td>0.0917</td><td>0.159</td></tr><tr><td rowspan="3">NES efficiency rate</td><td>Λ<sub>ml</sub><sup>th</sup></td><td>3.47</td><td>27,26</td><td>136.99</td></tr><tr><td>Λ<sub>ml,SF</sub><sup>num</sup></td><td>3.47</td><td>25,75</td><td>121.92</td></tr><tr><td>Λ<sub>ml,RM</sub><sup>num</sup></td><td>3.47</td><td>25,62</td><td>117.81</td></tr><tr><td>Relative error</td><td>100 × <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">Λ<sub>ml</sub><sup>th</sup> - Λ<sub>ml,RM</sub><sup>num</sup></div><div style="padding: 0 5px;">Λ<sub>ml,RM</sub><sup>num</sup></div></div></td><td>4.2 %</td><td>6.4 %</td><td>16.28 %</td></tr></tbody></table>

<figure><img src="image_69.png" /><figcaption>**FIGURE A.6.** (a) Sketch of the two DOFs airfoil (primary system) coupled to three one-DOF NESs. (b) Zoom on the *m*-th NES (with *m* = 1, ..., *M*).</figcaption></figure>

where, similarly as in Sect. A.4.1, Eqs. (A.72) correspond to Eqs. (D.3) and the corresponding Eqs. (D.5) and Eqs. (D.10) hereafter referred indistinctly as the Reference Model can be easily deduced as well as the slow-flow model corresponding to Eqs. (A.31).

The bifurcation parameter under consideration is now the reduced speed of the flow $\Theta$.

In the first instance, a model with three NESs ($M=3$) is considered and we use the following parameters of the primary system: $\xi_x = 1$, $\xi_\varphi = 1$, $\xi_x = 0.01$, $\xi_\varphi = 0.01$, $s_\varphi = 0.2$, $r_\varphi = 0.5$, $\Omega = 0.5$, $\eta = 0$, $\Lambda = 2\pi$ and $\mu = 1/(10\pi)$. Two sets of NESs parameters are used, hereafter referred as configuration 1 and configuration 2. The configuration 1 is defined with $\gamma_{h_1} = 0.4$, $\gamma_{h_2} = 0.3$, $\gamma_{h_3} = 0.2$, $\alpha_{h_1} = 7$, $\alpha_{h_2} = 6$, $\alpha_{h_3} = 6$, $\delta_1 = -0.9$, $\delta_2 = -0.8$ and $\delta_3 = -0.7$ which considers almost identical NESs. The configuration 2 is defined with $\gamma_{h_1} = 0.35$, $\gamma_{h_2} = 0.3$, $\gamma_{h_3} = 0.25$, $\alpha_{h_1} = 7$, $\alpha_{h_2} = 6$, $\alpha_{h_3} = 1$, $\delta_1 = -0.9$, $\delta_2 = -0.85$ and $\delta_3 = -0.8$ where the third NES has a small cubic stiffness compared to those of other NESs. Moreover, for both configurations $m_{h_1} = m_{h_2} = m_{h_3} = 1$ and several values for $\epsilon$ will be considered.

With this set of numerical values, the airfoil model (i.e. the primary system) undergoes an Hopf bifurcation at $\Theta_{Hopf}^{wo} = 0.933$ due to a mode coalescence phenomenon. If $\Theta < \Theta_{Hopf}^{wo}$ the system is stable and if $\Theta > \Theta_{Hopf}^{wo}$ the system is unstable. As previously, in unstable situations, the system reaches an harmful LCO we want to mitigate using the three NESs. We focus again on the possible mitigation regimes which can occur when the primary system is unstable.

First, in the Fig. A.7 the CM $S$, computed from (A.47) and (A.52), is depicted in the $(r_1, r_2, r_3)$-space with their stable parts (in gray), i.e. the $S_k$ (with $k = 1, ..., 2^M = 8$), and unstable parts (in black). The
stable and unstable fixed points, computed from (B.43) and (A.69), are again represented by magenta ans blue points respectively. We can see in each of the four figures a stable fixed point with a large value on the unbounded straight line (i.e. $S_{2M} = S_8$). As in Fig. A.4, this fixed point corresponds to a large amplitude LCO and if it is reached, we are in a harmful situations.

The theoretical trajectory of the slow-flow is determined for the four values of the bifurcation parameters $\Theta$ (in green) and compared to the trajectory obtained from the numerical simulation of the slow-flow (in red). The Fig. A.7(a) shows a harmless situation through a periodic regime, i.e. a stable fixed point of the slow-flow (in magenta) which is not in $S_8$ (here it is in $S_1$) is reached. In Figs. A.7(b) and A.7(c) harmless situations through SMRs are shown. Indeed, relaxation oscillations of the slow-flow are observed with respect to $r_1$ in Fig. A.7(b) and to $r_1$ and $r_2$ in Fig. A.7(c). Finally, because a stable fixed point is reached in $S_8$, the Fig. A.7(d) corresponds to an harmful situation.

Similar observations as in Fig. A.4 can be expressed: (1) during the slow epochs on $S$ the theoretical and numerical trajectories of the slow-flow are again almost superimposed and (2) due to the 0-order approximation of the asymptotic approach, differences are observed during fast epochs. Indeed, in the latter case, unlike the theoretical trajectory, the numerical trajectory is not at each jump parallel to only one direction of the $(r_1, r_2, r_3)$-space. Again, figures similar to Fig. A.7 using $\epsilon = 10^{-4}$ would show a perfect agreement between theoretical and numerical trajectories.

The Fig. A.8 shows the direct numerical integration of the Reference Model (blue lines) and of the slow-flow (red lines) for situations corresponding to Figs. A.7(c) and A.7(d). In Fig. A.8(a), $w_1(t)$ and $w_2(t)$ are strongly amplitude modulated signals (i.e. SMRs) which corresponds to relaxation oscillations of the slow-flow in both $r_1$ and $r_2$ directions. A large (harmful) amplitude periodic regime corresponding to a large stable fixed point of the slow-flow is observed in Fig. A.8(b).

Now the theoretical value of the mitigation limit $\Theta_{ml}^{th}$ is compared to the numerical estimation of the mitigation limit "measured" on the graphs of the maximum steady-state amplitude obtained from numerical simulations of the slow-flow $\Theta_{ml,SF}^{num}$ and from numerical simulations of the Reference model $\Theta_{ml,RM}^{num}$.

Again, the comparison is also performed in term of NES efficiency rate $\Lambda$ (in %), now defined with respect to the parameter $\Theta$

$$ \Lambda_{ml}^{th} = 100 \times \frac{\Theta_{ml}^{th} - \Theta_{Hopf}^{wo}}{\Theta_{Hopf}^{wo}}, \qquad (A.73a) $$

$$ \Lambda_{ml,SF}^{num} = 100 \times \frac{\Theta_{ml,SF}^{num} - \Theta_{Hopf}^{wo}}{\Theta_{Hopf}^{wo}}, \qquad (A.73b) $$

$$ \Lambda_{ml,RM}^{num} = 100 \times \frac{\Theta_{ml,RM}^{num} - \Theta_{Hopf}^{wo}}{\Theta_{Hopf}^{wo}}. \qquad (A.73c) $$

The comparison is presented in the Tab. A.2 for both configurations 1 and 2. Moreover, numerical bifurcation diagrams for the variables $v_1$ and $v_2$ (calculated from (A.10)) of the reference model are shown in Fig. A.9 (in red). The diagrams are compared to diagrams obtained for the primary system alone (in blue). Each of them is obtained from direct numerical integration of the corresponding system of equations and plotting the maximum steady-state amplitude, denoted $A_{v_n}$ ($n = 1, 2$), of each variable as functions of $\Theta$.

Once more, similar conclusions as in Sect. A.4.1 can be made from Tab. A.2. Indeed, Tab. A.2(a) shows that the mitigation limits obtained from numerical simulations of the slow-flow and of the reference model are very close and a difference of the order of magnitude of $\epsilon$ is observed between theoretical and numerical values. In term of NES efficiency rate (see Tab. A.2(b)), we can see that in efficient situations (i.e. $\epsilon = 10^{-2}$) the theoretical prediction can lead to a relatively large overestimation of the NES efficiency rate, for the same reasons as above.

As a final step, the influence of the number of NESs is examined. The results are shown in Tab. A.3 in which the configuration 1 is considered for a number of NESs from $M = 1$ to $M = 4$. In the case of $M = 4$ a fourth NES is added with the following parameters: $\gamma_{h_4} = 0.1$, $\alpha_{h_4} = 6$ and $\delta_4 = 0.6$. We can see that the prediction is relatively robust with respect to number of NESs even if a global slight increase of the relative error between theoretical and numerical efficiency rates is observed when $M$ increases.

## A.5 CONCLUSION

In this paper, dynamic instability mitigation of a multi-DOF system using Nonlinear Energy Sinks (NESs) is considered. The steady-state regimes of the coupled system is classified into two categories depending on whether the dynamic instability is mitigated or not and therefore separating harmless situations from harmful
<figure><img src="image_70.png" /><figcaption>**FIGURE A.7.** Critical manifold in the $(r_1, r_2, r_3)$-space defined through Eq. (E.32) together with theoretical and numerical trajectories of the slow-flow (A.31). Configuration 1 is considered with $\epsilon = 10^{-3}$.</figcaption></figure>

situations. The value of a chosen bifurcation parameter which separates harmless situations from harmful situations is called mitigation limit.

The prediction method of steady-state regimes has three main steps. The first step is the diagonalization of the primary system using the so-called biorthogonal transformation. Assuming a primary system with only one unstable mode the diagonalized system is reduced ignoring the stable modes and keeping only the unstable mode. The second step consists in applying the complexification-averaging method for the purposes of obtaining the slow-flow of the reduced system. It appears that the slow-flow is a slow-fast system. In light of that fact, in the third step, performing an asymptotic analysis using geometric singular perturbation theory of slow-fast systems, the critical manifold and the fixed points of the slow-flow are computed. The analysis shows that the critical manifold of the system can be reduced to a one dimensional parametric curve evolving
<figure><img src="image_71.png" /><figcaption>**FIGURE A.8.** Direct numerical integration of the Reference Model (solid blue line) and of the slow-flow (red line). Parameters: the same as those used for (a) Fig. A.7(c) and (b) Fig. A.7(d).</figcaption></figure>

in a multidimensional space. The response regimes of the system (and therefore the mitigation limit of the NES setup) is predicted by locating attracting parts, repelling parts and fold points of the critical manifold together with position and stability properties of the fixed points of the slow-flow. It is worth noting that the number of NESs does not appear as a theoretical restriction of the prediction method.

Finally, the method is applied to the prediction of the mitigation limit of two self-excited mechanical systems coupled to a NES setup. The first one is the well-known Hultén's model which reproduces the typical dynamic behavior of friction systems (it can undergo mode coupling instabilities) and coupled to two ungrounded NESs. The second one is an airfoil model undergoing an aeroelastic instability and coupled to a multi-NES setup (up to four). Theoretical results are compared, for validation purposes, to direct numerical integration of the system. The comparison shows a good agreement and highlights the limits of the asymptotic approach.

As a perspective, a generalization of this analysis could be suitable in the case of a primary system undergoing a multi-instability, as can be seen in some finite element models of braking systems. Indeed, during the first step of the method, not only one but several unstable modes would be kept. However, as a result of the increase of the critical manifold dimension, the prediction of the steady-state regimes would require a more advanced mathematical procedure.
<table><caption>**TABLEAU A.2.** Comparison between numerical and theoretical estimations of (a) the mitigation limit and (b) the NES efficiency rate for several values of the little parameter $\epsilon$.</caption><thead><tr><th colspan="2"></th><th colspan="3">(a)</th></tr><tr><th colspan="2"></th><th colspan="3">Mitigation limit</th></tr><tr><th colspan="2"></th><th colspan="3">$\epsilon$</th></tr><tr><th colspan="2"></th><th colspan="3">$10^{-4}$</th></tr></thead><tbody><tr><th rowspan="3">Configuration 1</th><th>$\Theta_{ml}^{th}$</th><td>0.93448</td><td>0.9471</td><td>1.057</td></tr><tr><th>$\Theta_{ml,SF}^{num}$</th><td>0.93446</td><td>0.9460</td><td>1.021</td></tr><tr><th>$\Theta_{ml,RM}^{num}$</th><td>0.93444</td><td>0.9459</td><td>1.020</td></tr><tr><th rowspan="3">Configuration 2</th><th>$\Theta_{ml}^{th}$</th><td>0.93430</td><td>0.9453</td><td>1.038</td></tr><tr><th>$\Theta_{ml,SF}^{num}$</th><td>0.93428</td><td>0.9444</td><td>1.013</td></tr><tr><th>$\Theta_{ml,RM}^{num}$</th><td>0.93426</td><td>0.9442</td><td>1.013</td></tr></tbody></table>

<table><caption>**TABLEAU A.3.** Comparison between numerical and theoretical estimations of the mitigation limit and the NES efficiency rate for several number of NESs attached on the airfoil for $\epsilon = 10^{-3}$. Configuration 1 is considered for a number of NESs from $M = 1$ to $M = 3$. In the case of $M = 4$, a fourth NES is added with the following parameters: $\gamma_{h_4} = 0.1$, $\alpha_{h_1} = 6$ and $\delta_3 = 0.7$.</caption><thead><tr><th colspan="2"></th><th colspan="3">(b)</th></tr><tr><th colspan="2"></th><th colspan="3">NES efficiency rate</th></tr><tr><th colspan="2"></th><th colspan="3">$\epsilon$</th></tr><tr><th colspan="2"></th><th colspan="3">$10^{-4}$</th></tr></thead><tbody><tr><th rowspan="3">Configuration 1</th><th>$\Lambda_{ml}^{th}$</th><td>0.16</td><td>1.51</td><td>13.29</td></tr><tr><th>$\Lambda_{ml,SF}^{num}$</th><td>0.16</td><td>1.39</td><td>9.43</td></tr><tr><th>$\Lambda_{ml,RM}^{num}$</th><td>0.16</td><td>1.38</td><td>9.32</td></tr><tr><th rowspan="3">Configuration 2</th><th>$\Lambda_{ml}^{th}$</th><td>0.14</td><td>1.32</td><td>11.25</td></tr><tr><th>$\Lambda_{ml,SF}^{num}$</th><td>0.14</td><td>1.22</td><td>8.57</td></tr><tr><th>$\Lambda_{ml,RM}^{num}$</th><td>0.14</td><td>1.2</td><td>8.57</td></tr></tbody></table>

<table><thead><tr><th colspan="2"></th><th colspan="4">Number of NESs attached (M)</th></tr><tr><th colspan="2"></th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><th rowspan="3">Mitigation limit</th><th>$\Theta_{ml}^{th}$</th><td>0.9392</td><td>0.9440</td><td>0.9471</td><td>0.9486</td></tr><tr><th>$\Theta_{ml,SF}^{num}$</th><td>0.9389</td><td>0.9432</td><td>0.9460</td><td>0.9469</td></tr><tr><th>$\Theta_{ml,RM}^{num}$</th><td>0.9388</td><td>0.9430</td><td>0.9459</td><td>0.9467</td></tr><tr><th rowspan="3">NES efficiency rate</th><th>$\Lambda_{ml}^{th}$</th><td>0.66</td><td>1.17</td><td>1.51</td><td>1.64</td></tr><tr><th>$\Lambda_{ml,SF}^{num}$</th><td>0.63</td><td>1.08</td><td>1.39</td><td>1.47</td></tr><tr><th>$\Lambda_{ml,RM}^{num}$</th><td>0.62</td><td>1.06</td><td>1.38</td><td>1.45</td></tr><tr><th>Relative error</th><td>100 × <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">$\Lambda_{ml}^{th} - \Lambda_{ml,RM}^{num}$</div><div style="padding: 0 5px;">$\Lambda_{ml,RM}^{num}$</div></div></td><td>6.45 %</td><td>10.38 %</td><td>9.42 %</td><td>13.1 %</td></tr></tbody></table>

## A.6 APPENDIX A. THE NONLINEAR HULTÈN'S MODEL WITH TWO NESs

The system is similar to the one used in [a14], it is composed of the nonlinear Hultèn's model [184] coupled to two ungrounded NESs with masses $\tilde{m}_{h_m}$, linear stiffnesses $k_{h_m}$, damping coefficients $\tilde{c}_{h_m}$ and cubic stiffnesses $k_{h_m}^{NL}$ ($m = 1, 2$, see Fig. F.1). The motion equations of the system are

$$ \begin{aligned} \ddot{x}_1 + c_{x_1} \dot{x}_1 + k_{x_1} \tilde{x}_1 - \mu k_{x_2} \tilde{x}_2 + \epsilon (k_{x_1}^{NL} \tilde{x}_1^3 - \mu k_{x_2}^{NL} \tilde{x}_2^3) + \\ \tilde{c}_{h_1} (\dot{x}_1 - \tilde{h}_1) + k_{h_1}^{NL} (\tilde{x}_1 - \tilde{h}_1)^3 = 0 \end{aligned} \quad (A.74a) $$ $$ \ddot{x}_2 + c_{x_2} \dot{x}_2 + k_{x_2} \tilde{x}_2 + \mu k_{x_1} \tilde{x}_1 + \epsilon (\mu k_{x_1}^{NL} \tilde{x}_1^3 + k_{x_2}^{NL} \tilde{x}_2^3) + $$
<figure><img src="image_72.png" /><figcaption>**FIGURE A.9.** Comparison between numerical bifurcation diagrams of the airfoil model with and without NESs with respect to the reduced speed of the flow $\Theta$. Configuration 1 is considered with (a) $\epsilon = 10^{-4}$, (b) $\epsilon = 10^{-3}$ and (c) $\epsilon = 10^{-2}$.</figcaption></figure>

$$ \tilde{c}_{h2}(\dot{x}_2 - \dot{h}_2) + k_{h2}^{NL}(\tilde{x}_2 - \tilde{h}_2)^3 = 0 \quad (A.74b) $$

$$ \bar{m}_{h_m} \ddot{h}_m - \tilde{c}_{h_m} (\dot{x}_m - \dot{h}_m) + k_{h_m}^{NL} (\tilde{x}_m - \tilde{h}_m)^3 = 0 \quad (m = 1, 2) \quad (A.74c) $$

where the time derivative is denoted " · ". In this model the tangential force $F_T$ due to friction contact is assumed to be proportional to the normal force $F_N$ as given by Coulomb's law: $F_T = \mu F_N$, where $\mu$ is the friction coefficient.

Introducing the following notation: $\omega_n = \sqrt{k_{x_n}/m}$, $\Omega = \omega_2/\omega_1$, $\zeta_{x_n} = c_{x_n}/(m\omega_1)$, $\xi_{x_n} = k_{x_n}^{NL}/(m\omega_1^2)$ (with $n = 1, 2$) for the primary system parameters and $\tilde{m}_{h_m} = \tilde{m}_{h_m}/(m\omega_1^2)$, $\tilde{\zeta}_{h_m} = \tilde{c}_{h_m}/(m\omega_1)$ and $\xi_{x_m} = k_{x_m}^{NL}/(m\omega_1^2)$ (with $m = 1, 2$) for the NESs parameters, Eqs. (A.74) become

$$ x_1'' + \zeta_{x_1} x_1' + x_1 - \mu \Omega^2 x_2 + \epsilon (\xi_{x_1} x_1^3 - \mu \xi_{x_2} x_2^3) + \tilde{m}_{h_1} h_1'' = 0 \quad (A.75a) $$

$$ x_2'' + \zeta_{x_2} x_1' + \Omega^2 x_2 + \mu x_1 + \epsilon (\mu \xi_{x_1} x_1^3 + \xi_{x_2} x_2^3) + \tilde{m}_{h_2} h_2'' = 0 \quad (A.75b) $$

$$ h_m'' - \gamma_{h_m} (x_m' - h_m') - \alpha_{h_m} (x_m - h_m)^3 = 0 \quad (m = 1, 2) \quad (A.75c) $$

The system of equations (A.75) has the general form of Eq. (D.2) with

$$ \tilde{\mathbf{x}} = \begin{pmatrix} \tilde{x}_1 \\ \tilde{x}_2 \end{pmatrix}, \quad \tilde{\mathbf{h}} = \begin{pmatrix} \tilde{h}_1 \\ \tilde{h}_2 \end{pmatrix}, \quad (A.76) $$

$$ \tilde{\mathbf{M}} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}, \quad \tilde{\mathbf{C}} = \begin{bmatrix} \zeta_{x_1} & 0 \\ 0 & \zeta_{x_1} \end{bmatrix}, \quad \tilde{\mathbf{K}} = \begin{bmatrix} 1 & -\mu\Omega^2 \\ \mu & \Omega^2 \end{bmatrix}, \quad (A.77a) $$

$$ \tilde{\mathbf{G}} = \begin{bmatrix} \xi_x & 0 \\ 0 & \xi_\varphi \end{bmatrix}, \quad \tilde{\mathbf{B}} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \quad \text{and} \quad \mathbf{T} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}. \quad (A.77b) $$

The nonlinear vector function contains cubic nonlinearities

$$ \mathbf{g}^{\text{NL}}(\tilde{\mathbf{x}}) = \begin{pmatrix} \xi_{x_1} x_1^3 - \mu \xi_{x_2} x_2^3 \\ \mu \xi_{x_1} x_1^3 + \xi_{x_2} x_2^3 \end{pmatrix}, \quad (A.78) $$

and the nonlinearities of the NESs have the required form, i.e. $\mathbf{f}^{\text{NL}}(\mathbf{T}\tilde{\mathbf{x}} - \tilde{\mathbf{h}})$, with $\mathbf{f}^{\text{NL}}$ given by (A.5).

According to the method proposed in Sect. D.2.1, the small parameter $\epsilon$ is introduced to take into account small inertia and damping terms of the NESs

$$ \tilde{m}_{h_m} = \epsilon m_{h_m}, \quad \tilde{\zeta}_{h_m} = \epsilon \zeta_{h_m} \quad (m = 1, 2), \quad (A.79) $$

and the variables $\tilde{x}_n$ ($n = 1, 2$) and $\tilde{h}_m$ ($m = 1, 2$) are rescaled through $\epsilon$ as

$$ x_n = \frac{\tilde{x}_n}{\sqrt{\epsilon}} \quad (n = 1, 2) \quad h_m = \frac{\tilde{h}_m}{\sqrt{\epsilon}} \quad (m = 1, 2). \quad (A.80) $$
Consequently, the equations of motion have the following final form

$$x_1'' + \zeta_{x_1} x_1' + x_1 - \mu \Omega^2 x_2 + \epsilon (\xi_{x_1} x_1^3 - \mu \xi_{x_2} x_2^3) + \epsilon m_{h_1} h_1'' = 0 \quad (A.81a)$$

$$x_2'' + \zeta_{x_2} x_1' + \Omega^2 x_2 + \mu x_1 + \epsilon (\mu \xi_{x_1} x_1^3 + \xi_{x_2} x_2^3) + \epsilon m_{h_2} h_2'' = 0 \quad (A.81b)$$

$$h_m'' - \gamma_{h_m} (x_m' - h_m') - \alpha_{h_m} (x_m - h_m)^3 = 0 \quad (m = 1, 2) \quad (A.81c)$$

with $\gamma_{h_m} = \zeta_{h_m}/m_{h_m}$ and $\alpha_{h_m} = \xi_{h_m}/m_{h_m}$ ($m = 1, 2$).

## A.7 APPENDIX B. THE AIRFOIL-M-NESs MODEL

The two degrees of freedom airfoil model coupled to several NESs (see Fig. D.3) and used in Sect. D.5 is described in this appendix.

Assuming small angle, we derive the equations of motion of the airfoil-M-NESs model as

$$m \ddot{z} + S_{\varphi} \ddot{\varphi} + K_z \ddot{z} + K_z^{NL} \ddot{z}^3 + C_z \dot{z} + q P_a \Lambda (\dot{\varphi} + \dot{z}/U) + \sum_{m=1}^{M} \tilde{C}_{y_m} (\dot{z} - d_m \dot{\varphi} - \dot{y}_m) + K_{y_m}^{NL} (\tilde{z} - d_m \tilde{\varphi} - \tilde{y}_m)^3 = 0 \quad (A.82a)$$

$$I_{\varphi} \ddot{\varphi} + S_{\varphi} \ddot{z} + K_{\varphi} \ddot{\varphi} + K_{\varphi}^{NL} \ddot{\varphi}^3 + C_{\varphi} \dot{\varphi} - q e P_a \Lambda (\dot{\varphi} + \dot{z}/U) - \sum_{m=1}^{M} d_m \tilde{C}_{y_m} (\dot{z} - d_m \dot{\varphi} - \dot{y}_m) + d_m K_{y_m}^{NL} (\tilde{z} - d_m \tilde{\varphi} - \tilde{y}_m)^3 = 0 \quad (A.82b)$$

$$\tilde{m}_{y_m} \ddot{y}_m + \tilde{C}_{y_m} (d_m \dot{\varphi} + \dot{y}_m - \dot{z}) + K_{y_m}^{NL} (d_m \tilde{\varphi} + \tilde{y}_m - z)^3 = 0 \quad (m = 1, \dots, M) \quad (A.82c)$$

where again the time derivative is denoted "·". $b = c/2$ is the semichord length. A is the aerodynamic center, B the elastic axis, C the center of gravity of the airfoil. $e$ is the location aerodynamic center A measured from B (positive ahead of B). $m$ and $I_{\varphi}$ are the mass of the airfoil and its mass moment of inertia with respect to B. $S_{\varphi} = mBC$ is the mass unbalance in the airfoil. $K_z$ and $K_{\varphi}$ are the linear heave and pitch stiffnesses respectively whereas $K_z^{NL}$ and $K_{\varphi}^{NL}$ are the cubic heave and pitch stiffnesses. $C_z$ and $C_{\varphi}$ are the heave and pitch damping coefficients. $U$ is the constant and uniform flow speed around the airfoil and $q = \frac{1}{2} \rho_{\infty} U^2$ is the dynamic pressure where $\rho_{\infty}$ is the density of the flow. $P_a$ is the planform of the airfoil, $\Lambda$ is the lift curve slope and $d_m$ ($m = 1, \dots, M$) are the offset attachments of the NESs to the airfoil. Finally, $\tilde{m}_{y_m}$, $\tilde{C}_{y_m}$ and $K_{y_m}^{NL}$ ($m = 1, \dots, M$) are the masses, the damping coefficients and the cubic stiffnesses of the NESs respectively.

For convenience, Eq. (D.66) is written in dimensionless form

$$\tilde{x}'' + s_{\varphi} \tilde{\varphi}'' + \Omega^2 \tilde{x} + \xi_x \tilde{x}^3 + \zeta_x \tilde{x}' + \mu \Lambda \Theta (\Theta \varphi + \tilde{x}') + \sum_{m=1}^{M} \tilde{m}_{h_m} \tilde{h}_m'' = 0 \quad (A.83a)$$

$$r_{\varphi}^2 \tilde{\varphi}'' + s_{\varphi} \tilde{x}'' + r_{\varphi}^2 \tilde{\varphi} + \xi_{\varphi} \tilde{\varphi}^3 + \zeta_{\varphi} \tilde{\varphi}' - \eta \mu \Lambda \Theta (\Theta \varphi + \tilde{x}') - \sum_{m=1}^{M} \delta_m \tilde{m}_{h_m} \tilde{h}_m'' = 0 \quad (A.83b)$$

$$\tilde{m}_{h_m} \tilde{h}_m'' + \tilde{\zeta}_{h_m} (\delta_m \tilde{\varphi}' + \tilde{h}_m' - \tilde{x}') + \xi_{h_m} (\delta_m \tilde{\varphi} + \tilde{h}_m - \tilde{x})^3 = 0 \quad (m = 1, \dots, M) \quad (A.83c)$$

where the time $t$ has been replaced by the dimensionless time $\tau = \omega_{\varphi} t$ (with $\omega_{\varphi} = \sqrt{K_{\varphi}/I_{\varphi}}$) and time derivative $\partial/\partial\tau$ is denoted "·". The dimensionless displacements are defined by $x = z/b$, $h = y/b$. Moreover, $\delta_m = d_m/b$ ($m = 1, \dots, M$), $s_{\varphi} = S_{\varphi}/(mb)$, $\Omega = \omega_z/\omega_{\varphi}$ (with $\omega_z = \sqrt{K_z/m}$) and $\eta = e/b$. $r_{\varphi} = \sqrt{I_{\varphi}/(mb^2)}$ is the radius of gyration of the cross section of the wing. The dimensionless nonlinear stiffnesses and damping coefficients are $\xi_x = K_z^{NL} b^2/(m\omega_{\varphi}^2)$, $\zeta_x = C_z/(m\omega_{\varphi})$, $\xi_{\varphi} = K_{\varphi}^{NL} / (mb^2\omega_{\varphi}^2)$, $\zeta_{\varphi} = C_{\varphi}/(mb^2\omega_{\varphi})$, $\xi_{h_m} = K_{y_m}^{NL} b^2/(m\omega_{\varphi}^2)$, $\tilde{\zeta}_{h_m} = \tilde{C}_{y_m}/(m\omega_{\varphi})$. The mass and density ratios are defined by $\tilde{m}_{h_m} = \tilde{m}_{y_m}/(m\omega_{\varphi}^2)$ and $\mu = bP_a\rho_{\infty}/(2m)$ respectively. Finally, the bifurcation parameter under consideration is the reduced speed of the flow $\Theta = U/(b\omega_{\varphi})$.

The system of equations (D.67) has the general form of Eq. (D.2) with

$$\tilde{\mathbf{x}} = \begin{pmatrix} \tilde{x} \\ \tilde{\varphi} \end{pmatrix}, \quad \tilde{\mathbf{h}} = \begin{pmatrix} \tilde{h}_1 \\ \tilde{h}_2 \\ \tilde{h}_3 \end{pmatrix}, \quad (A.84)$$
$$ \tilde{\mathbf{M}} = \begin{bmatrix} 1 & s_{\varphi} \\ s_{\varphi} & r_{\varphi}^2 \end{bmatrix}, $$

$$ \tilde{\mathbf{C}} = \begin{bmatrix} \zeta_x + \mu\Theta & 0 \\ -\eta\mu\Theta & \zeta_{\varphi} \end{bmatrix}, $$

$$ \tilde{\mathbf{K}} = \begin{bmatrix} \Omega^2 & \mu\Theta^2 \\ 0 & r_{\varphi}^2 - \eta\mu\Theta^2 \end{bmatrix}, \qquad (\text{A.85a}) $$

$$ \tilde{\mathbf{G}} = \begin{bmatrix} \xi_x & 0 \\ 0 & \xi_{\varphi} \end{bmatrix}, $$

$$ \tilde{\mathbf{B}} = \begin{bmatrix} 1 & \cdots & 1 \\ -\delta_1 & \cdots & -\delta_M \end{bmatrix} \quad \text{and} $$

$$ \mathbf{T} = \begin{bmatrix} 1 & -\delta_1 \\ \vdots & \vdots \\ 1 & -\delta_M \end{bmatrix}. \qquad (\text{A.85b}) $$

The expression of the nonlinear vector function $\mathbf{g}^{\text{NL}}$ is given by

$$ \mathbf{g}^{\text{NL}}(\tilde{\mathbf{x}}) = \begin{pmatrix} \xi_x \tilde{x}^3 \\ \xi_{\varphi} \tilde{\varphi}^3 \end{pmatrix}, \qquad (\text{A.86}) $$

and again $\mathbf{f}^{\text{NL}}(\mathbf{T}\tilde{\mathbf{x}} - \tilde{\mathbf{h}})$, with $\mathbf{f}^{\text{NL}}$ given by (A.5).

Following Eqs. (A.6) and (A.7), we define

$$ \tilde{m}_{h_m} = \epsilon m_{h_m}, \quad \tilde{\zeta}_{h_m} = \epsilon \zeta_{h_m} \quad (m = 1, \dots, M) \qquad (\text{A.87}) $$

and

$$ x = \frac{\tilde{x}}{\sqrt{\epsilon}}, \quad \varphi = \frac{\tilde{\varphi}}{\sqrt{\epsilon}}, \quad h_m = \frac{\tilde{h}_m}{\sqrt{\epsilon}} \quad (m = 1, \dots, M) \qquad (\text{A.88}) $$

to obtain the final form of the equations of motion

$$ x'' + s_{\varphi}\varphi'' + \Omega^2 x + \epsilon\xi_x x^3 + \zeta_x x' + \mu\Lambda\Theta(\Theta\varphi + x') + \epsilon \sum_{m=1}^{M} m_{h_m} h_m'' = 0 \qquad (\text{A.89a}) $$

$$ r_{\varphi}^2 \varphi'' + s_{\varphi} x'' + r_{\varphi}^2 \varphi + \epsilon \xi_{\varphi} \varphi^3 + \zeta_{\varphi} \varphi' - \eta \mu \Lambda \Theta (\Theta \varphi + x') - \epsilon \sum_{m=1}^{M} \delta_m m_{h_m} h_m'' = 0 \qquad (\text{A.89b}) $$

$$ h_m'' + \gamma_{h_m} (\delta_m \varphi' + h_m' - x') + \alpha_{h_m} (\delta_m \varphi + h_m - x)^3 = 0 \quad (m = 1, \dots, M) \qquad (\text{A.89c}) $$

with $\gamma_{h_m} = \zeta_{h_m}/m_{h_m}$ and $\alpha_{h_m} = \xi_{h_m}/m_{h_m}$ ($m = 1, \dots, M$).
# Annexe B

Communications in Nonlinear Science and Numerical Simulation 95 (2021) : 105623 [a9]

## DYNAMIC BEHAVIOR ANALYSIS OF A MECHANICAL SYSTEM WITH TWO UNSTABLE MODES COUPLED TO A SINGLE NONLINEAR ENERGY SINK

Baptiste BERGEOT, Sergio BELLIZZI and Sébastien BERGER

### ABSTRACT

This paper investigates a problem of passive mitigation of vibratory instabilities caused by two unstable modes by means of a single nonlinear energy sink (NES). For this purpose, a linear four-degree-of-freedom (DOF) primary structure having two unstable modes (reproducing the typical dynamic behavior of a friction system) and undergoing, as it is linear, unbounded motions when it is unstable, is coupled to a NES. In this work, the NES involves an essentially cubic restoring force and a linear damping force. We are interested in characterizing analytically the response regimes resulting from the coupling of the two unstable linear modes of the primary structure and the nonlinear mode of the NES. To this end, from a suitable rescaling of the governing equations of the coupled system in which the dynamics of the primary structure is reduced to its unstable modal coordinates, the complexification-averaging method is applied. The resulting averaged system appears to be a fast-slow system with four fast variables and two slow ones related to the two unstable modes of the primary structure. The critical manifold of the averaged dynamics is obtained through the geometric singular perturbation theory and appears as a two-dimensional parametric surface (with respect to two of the four fast variables) which evolves in the whole six-dimensional variable space. The asymptotic analysis reveals that the NES attachment can produce some bounded responses and suggests that the system may have simultaneous stable attractors. Numerical simulations complement the study, highlighting a possible competition between stable attractors and allowing us to investigate their basins of attraction. In each considered situation, a good agreement has been observed between theoretical results and numerical simulations, which validates the proposed asymptotic analysis.

**Keywords** : Nonlinear energy sink, Multi-instabilities, Relaxation oscillations, Multiple-scale analysis

### B.1 INTRODUCTION

A nonlinear energy sink (NES) is a vibration absorber used for passive shock and vibration mitigation of undesirable oscillations of a primary structure induced by external, parametric or self-excitations. A 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 and a viscous linear damper. Properly designed, a NES can be tuned to any frequency to perform an energy transfer from the primary structure towards itself. This phenomenon, called targeted energy transfer (TET), is realized through nonlinear mode bifurcations and it can be described as resonance captures. Reviews of these concepts can be found in [206] and in a part of [129].

The use of NES to control or suppress dynamic instabilities has been widely studied. The first work reported in this field concerns the suppression of limit cycle oscillations of the Van der Pol oscillator [122]. Bifurcations of self-excitation regimes in a Van der Pol (respectively Van der Pol-Duffing) oscillator coupled with a NES have been considered in [75] (respectively [59]). Numerous papers [78, 94, 121, 123, 124, 131, 156, 207] have been dedicated to the suppression of aeroelastic instabilities. The mitigation, by means of a NES, of vortex-induced vibrations resulting from the nonlinear interaction of a laminar flow and a rigid circular cylinder has been described in [47, 201]. The passive suppression of helicopter ground resonance instabilities has been considered in [a10, a12]. Chatter control in machine tool vibrations has been studied in [82, 142]. Mitigation of friction-induced vibrations due to mode coupling instabilities in a friction system has been investigated in [a14]. A network of parallel NES coupled to a Van der Pol oscillator has been considered in [a7]. The prediction of the steady-state regimes of a multi-degree-of-freedom dynamical system having one unstable mode coupled to a set of NES has been performed in [a8].
The presence of a NES modifies the position of the bifurcation point and the amplitude of the bifurcated branches. It has been shown that NES can mitigate limit cycle oscillations (LCO). Several possible scenarios have been identified including complete suppression, mitigation through periodic response and mitigation through strongly modulated response (SMR). The results have been obtained under the main assumption that the primary system has only one unstable mode. This assumption is the common thread in all previously referenced works.

Mitigations involving two simultaneous resonances have been considered in the cases of external excitations or transient regimes. TET by means of a multi-resonance capture in a 2-DOF system under two different harmonic excitations has been investigated numerically [13, 159, 160] and experimentally [45]. It has been shown that two resonances of the system are possible simultaneously resulting in vibration reduction around the two resonance frequencies. Regarding transient responses, multi-frequency TET cascades have been studied in [113, 202, 204]. Moreover, it has been shown in [68, 95] that a bistable NES is able to efficiently mitigate oscillations of the primary system for more than one vibration mode.

In this paper, we consider the case of a primary system having two unstable modes and we investigate the mitigation of the two instabilities by means of a single NES. To the knowledge of the authors, this case has never been considered in the past. As for the primary system, the retained model, a simple 4-DOF friction system involving two unstable modes, has been introduced in [53] and used in [52] to illustrate different stability behaviors characterized by mode coupling having coalescence patterns. Based on the procedure presented in [a8], we provide an analytic framework for characterizing the response regimes.

The paper is organized as follows. In Sect. B.2, the system under study is presented. It consists in two 2-DOF Hulten systems connected by one spring and one damper resulting in a 4-DOF system on which one ungrounded NES is fixed along one of the DOF. In Sect. B.3, applying a complexification-averaging method, the analyses of the critical manifold, fixed points and relaxation oscillations regimes are addressed. Some numerical examples are given and discussed in order to compare and to complete the asymptotic analysis in Sect. B.4. Finally, the paper is concluded in Sect. F.8.

## B.2 THE SYSTEM UNDER STUDY

### B.2.1 The mechanical model

The system considered here is shown in Fig. B.1. It is composed of a primary system coupled to a NES in an ungrounded configuration. The primary system is defined by two 2-DOF Hulten systems with respective masses $m_1$ and $m_2$ which are linearly coupled by a spring of stiffness $k_a$, and a viscous damper of damping coefficient $c_a$. The two masses are arranged against three moving belts. It is assumed that the belts run at a constant speed, that the masses and the belts are always in contact and that the relative velocities between the masses and the belts are always positive so that the directions of the tangential fiction forces do not change. Hence, the tangential forces $F_T$ due to friction contacts are proportional to the normal forces $F_N$ as given by the Coulomb's law : $F_T = \gamma F_N$, where $\gamma$ is the friction coefficient. This 4-DOF linear model was introduced in [53] and used in [52] to illustrate different stability behaviors characterized by mode coupling having coalescence patterns. The number of DOF of the primary system corresponds to the minimal number of DOF which can be used to investigate the case of multi-instabilities.

The NES is simply modeled as a 1-DOF nonlinear oscillator characterized by a mass $m_h$, no linear stiffness, a cubic stiffness $k_h^{NL}$ and a damping coefficient $c_h$.

The system is described by the following equations of motion

$$m_1 \frac{d^2 \tilde{x}_1}{d\tilde{\tau}^2} + c_1 \frac{d\tilde{x}_1}{d\tilde{\tau}} + c_a \left( \frac{d\tilde{x}_1}{d\tilde{\tau}} - \frac{d\tilde{x}_3}{d\tilde{\tau}} \right) + k_1 \tilde{x}_1 + k_a (\tilde{x}_1 - \tilde{x}_3) + \gamma k_2 \tilde{x}_2 + m_h \frac{d^2 \tilde{h}}{d\tilde{\tau}^2} = 0 \quad (B.1a)$$

$$m_1 \frac{d^2 \tilde{x}_2}{d\tilde{\tau}^2} + c_2 \frac{d\tilde{x}_2}{d\tilde{\tau}} + k_2 \tilde{x}_2 - \gamma k_1 \tilde{x}_1 = 0 \quad (B.1b)$$

$$m_2 \frac{d^2 \tilde{x}_3}{d\tilde{\tau}^2} + c_3 \frac{d\tilde{x}_3}{d\tilde{\tau}} + c_a \left( \frac{d\tilde{x}_3}{d\tilde{\tau}} - \frac{d\tilde{x}_1}{d\tilde{\tau}} \right) + k_3 \tilde{x}_3 + k_a (\tilde{x}_3 - \tilde{x}_1) + \gamma k_4 \tilde{x}_4 = 0 \quad (B.1c)$$

$$m_2 \frac{d^2 \tilde{x}_4}{d\tilde{\tau}^2} + c_4 \frac{d\tilde{x}_4}{d\tilde{\tau}} + k_4 \tilde{x}_4 - \gamma k_3 \tilde{x}_3 = 0 \quad (B.1d)$$
<figure><img src="image_73.png" /><figcaption>FIGURE B.1. The phenomenological system under study.</figcaption></figure>

$$m_h \frac{d^2 \tilde{h}}{d \tilde{\tau}^2} + c_h \left( \frac{d \tilde{h}}{d \tilde{\tau}} - \frac{d \tilde{x}_1}{d \tilde{\tau}} \right) + k_h^{NL} (\tilde{h} - \tilde{x}_1)^3 = 0 \quad (B.1e)$$

where $\tilde{x}_1$ and $\tilde{x}_2$ (respectively $\tilde{x}_3$ and $\tilde{x}_4$) are the displacements of the primary mass $m_1$ (respectively $m_2$) and $\tilde{h}$ is the displacement of the NES.

Introducing the following notations

$$\omega_1 = \sqrt{\frac{k_1}{m_1}}, \omega_2 = \sqrt{\frac{k_2}{m_1}}, \omega_3 = \sqrt{\frac{k_3}{m_2}}, \omega_4 = \sqrt{\frac{k_4}{m_2}},$$

$$\eta_i = \frac{c_1}{\omega_1 m_1}, \eta_2 = \frac{c_2}{\omega_2 m_1}, \eta_3 = \frac{c_3}{\omega_3 m_2}, \eta_4 = \frac{c_4}{\omega_4 m_2},$$

$$\Omega_i = \frac{\omega_i}{\omega_1} \quad \text{for } i = 2, \dots, 4,$$

$$\sigma = \frac{m_2}{m_1}, \Omega_a = \frac{k_a}{\omega_1 m_1}, \eta_a = \frac{c_a}{\omega_1 m_1},$$

$$\tilde{\mu} = \frac{c_h}{\omega_1 m_1}, \alpha = \frac{k_h^{NL}}{\omega_1^2 m_1}, \epsilon = \frac{m_h}{m_1},$$

we can then rescale the variables $\tilde{x}_i$, for $i = 1$ to 4 and $\tilde{h}$ as

$$x_i = \frac{\tilde{x}_i}{\sqrt{\epsilon}} \quad \text{for } i = 1, \dots, 4 \quad \text{and} \quad h = \frac{\tilde{h}}{\sqrt{\epsilon}}, \quad (B.2)$$

and change the time scale from $\tilde{\tau}$ to $t = \omega_1 \tilde{\tau}$. Denoting as "." the derivative with respect to time $t$, the system of equations (B.1) takes the following final form

$$\ddot{x}_1 + (\eta_1 + \eta_a)\dot{x}_1 - \eta_a\dot{x}_3 + (1 + \Omega_a^2)x_1 + \gamma\Omega_2^2 x_2 - \Omega_a^2 x_3 + \epsilon\ddot{h} = 0 \quad (B.3a)$$

$$\ddot{x}_2 + \eta_2 \Omega_2 \dot{x}_2 - \gamma x_1 + \Omega_2^2 x_2 = 0 \quad (B.3b)$$

$$\ddot{x}_3 - \frac{\eta_a}{\sigma}\dot{x}_1 + \left(\eta_3\Omega_3 + \frac{\eta_a}{\sigma}\right)\dot{x}_3 - \frac{\Omega_a^2}{\sigma}x_1 + \left(\Omega_3^2 + \frac{\Omega_a^2}{\sigma}\right)x_3 + \gamma\Omega_4^2 x_4 = 0 \quad (B.3c)$$

$$\ddot{x}_4 + \eta_4 \Omega_4 \dot{x}_4 - \gamma \Omega_3^2 x_3 + \Omega_4^2 x_4 = 0 \quad (B.3d)$$

$$\ddot{h} + \mu (\dot{h} - \dot{x}_1) + \alpha (h - x_1)^3 = 0 \quad (B.3e)$$
where $\tilde{\mu}$ has been rescaled as $\tilde{\mu} = \epsilon\mu$.

From now on we assume that $\epsilon$ is small. As it also implies that $\tilde{\mu}$ is small, this assumption is satisfied if the mass ratio and the damping coefficient of the NES are small.

### B.2.2 The reduced model

To simplify the following calculations, it is convenient to introduce new coordinates in (B.3) as

$$u_1 = x_1 + \epsilon h, \quad u_i = x_i \quad \text{for } i = 2, \dots, 4 \quad \text{and} \quad v = x_1 - h \qquad (B.4)$$

leading, after a first-order Taylor expansion around $\epsilon = 0$, to the following equations

$$\ddot{u}_1 + (\eta_1 + \eta_a)\dot{u}_1 - \eta_a\dot{u}_3 + (1 + \Omega_a^2)u_1 + \gamma\Omega_2^2 u_2 - \Omega_a^2 u_3 + \epsilon\left[-(\eta_1 + \eta_a)\dot{u}_1 + (\eta_1 + \eta_a)\dot{v} - (1 + \Omega_a^2)u_1 + (1 + \Omega_a^2)v\right] = 0 \qquad (B.5a)$$

$$\ddot{u}_2 + \eta_2 \Omega_2 \dot{u}_2 - \gamma u_1 + \Omega_2^2 u_2 + \epsilon \gamma (u_1 - v) = 0 \qquad (B.5b)$$

$$\ddot{u}_3 - \frac{\eta_a}{\sigma}\dot{u}_1 + \left(\eta_3\Omega_3 + \frac{\eta_a}{\sigma}\right)\dot{u}_3 - \frac{\Omega_a^2}{\sigma}u_1 + \gamma\Omega_4^2 u_4 + \left(\Omega_3^2 + \frac{\Omega_a^2}{\sigma}\right)u_3 + \frac{\epsilon}{\sigma}\left[\eta_a\dot{u}_1 - \eta_a\dot{v} + \Omega_a^2 u_1 - \Omega_a^2 v\right] = 0 \qquad (B.5c)$$

$$\ddot{u}_4 + \eta_4 \Omega_4 \dot{u}_4 - \gamma \Omega_3^2 u_3 + \Omega_4^2 u_4 = 0 \qquad (B.5d)$$

$$\ddot{v} + (\eta_1 + \eta_a)\dot{u}_1 - \eta_a\dot{u}_3 + \mu\dot{v} + (1 + \Omega_a^2)u_1 + \gamma\Omega_2^2 u_2 - \Omega_a^2 u_3 + \alpha v^3 + \epsilon\left[-(\eta_1 + \eta_a)\dot{u}_1 + (\eta_1 + \eta_a)\dot{v} + \mu\dot{v} - (1 + \Omega_a^2)u_1 + (1 + \Omega_a^2)v + \alpha v^3\right] = 0. \qquad (B.5e)$$

Then, to capture the essential features of the double instability of the primary system, a biorthogonal transformation is applied to the variables characterizing the primary system after writing Eqs. (B.5a) to (B.5d) in state-space form:

$$\dot{\mathbf{y}} = \mathbf{A}\mathbf{y} + \epsilon (\mathbf{D}_1\mathbf{y} + \mathbf{D}_2v + \mathbf{D}_3\dot{v}) \qquad (B.6)$$

where $\mathbf{y} = (u_1, \dots, u_4, \dot{u}_1, \dots, \dot{u}_4)^T$. The matrices $\mathbf{A}$ and $\mathbf{D}_1$ and the vectors $\mathbf{D}_2$ and $\mathbf{D}_3$ are not explicitly given but can be easily obtained from (B.5a-B.5d). The non symmetric matrix $\mathbf{A}$ characterizes the dynamics of the primary system.

Let us consider the right and left eigenvector matrices

$$\mathbf{R} = [\mathbf{r}_1 \ \mathbf{r}_1^* \ \dots \ \mathbf{r}_4 \ \mathbf{r}_4^*] \quad \text{and} \quad \mathbf{L} = [\mathbf{l}_1 \ \mathbf{l}_1^* \ \dots \ \mathbf{l}_2 \ \mathbf{l}_4^*] \qquad (B.7)$$

corresponding to the following right and left eigenvalue problems with respect to the matrix $\mathbf{A}$

$$\mathbf{A}\mathbf{R} = \mathbf{R}\mathbf{\Lambda} \quad \text{and} \quad \mathbf{A}^T\mathbf{L} = \mathbf{L}\mathbf{\Lambda}. \qquad (B.8)$$

where $\mathbf{\Lambda} = \text{diag}(\lambda_1, \lambda_1^*, \dots, \lambda_4, \lambda_4^*)$ denotes the diagonal matrix of the eigenvalues. The superscript $(\cdot)^*$ denotes the complex conjugate.

From now on we assume that the primary system is unstable with two pairs of conjugate eigenvalues having small positive real parts, we arrange the eigenvalues so that $\lambda_1^R > 0$ and $\lambda_2^R > 0$ where $\lambda_i = \lambda_i^R + j\lambda_i^I$ for $i = 1, 2$ and we assume that

$$\lambda_1^R = \epsilon \rho_1 \quad \text{and} \quad \lambda_2^R = \epsilon \rho_2 \qquad (B.9)$$

with $\rho_1$ and $\rho_2$ of the order of $O(1)$. The two other pairs of eigenvalues have negative real parts.

We introduce the biorthogonal transformation

$$\mathbf{y} = \mathbf{R}\mathbf{q} \quad \text{with} \quad \mathbf{q} = \mathbf{L}^T\mathbf{y} \qquad (B.10)$$

where the vector $\mathbf{q} = (q_1, q_1^*, \dots, q_4, q_4^*)^T$ is constituted of four pairs of complex conjugates. The matrices $\mathbf{R}$ and $\mathbf{L}$ have been chosen such that

$$\mathbf{L}^T\mathbf{R} = \mathbf{I} \quad \text{and} \quad \mathbf{L}^T\mathbf{A}\mathbf{R} = \mathbf{\Lambda}. \qquad (B.11)$$

where $\mathbf{I}$ denotes the identity matrix.

We can deduce that Eqs. (B.5) take the following formal form with respect to the new variables $\mathbf{q}$:

$$\dot{\mathbf{q}} - \mathbf{\Lambda}\mathbf{q} - \epsilon\mathbf{f}(\mathbf{q}, v, \dot{v}) = 0 \qquad (B.12a)$$

$$\ddot{v} + g(v, \dot{v}, \mathbf{q}, \epsilon) = 0 \qquad (B.12b)$$
where $\mathbf{f} = \mathbf{L}^T (\mathbf{D}_1 \mathbf{R}\mathbf{q} + \mathbf{D}_2 v + \mathbf{D}_3 \dot{v})$. Again, the function $g$ is not explicitly given but can be easily obtained from (B.5e).

Neglecting the coupling term $\epsilon \mathbf{f}(\mathbf{q}, v, \dot{v})$, Eq. (B.12a) reduces to $\dot{\mathbf{q}} - \mathbf{\Lambda}\mathbf{q} = 0$ and the variables $q_3, q_3^*, q_4$ and $q_4^*$, related to the stable modes, vanish (i.e. tend to zero). When the coupling term $\epsilon \mathbf{f}(\mathbf{q}, v, \dot{v})$, which is of the order of magnitude of $\epsilon$, is taken into account, the variables $q_3, q_3^*, q_4$ and $q_4^*$ do not vanish but their contributions remain minimal and can be neglected [a8, 78]. Consequently, all the terms related to $q_3, q_3^*, q_4$ and $q_4^*$ will be omitted in the following developments.

Eq. (B.12) may therefore be simplified, retaining only the components $q_1, q_1^*, q_2$ and $q_2^*$ associated with the unstable eigenvalues. The resulting equations take the (reduced) form

$$ \epsilon \mathbf{L}_1^T [\mathbf{D}_1 (\mathbf{R}_1 \cdot q_1 + \mathbf{R}_2 \cdot q_1^* + \mathbf{R}_3 \cdot q_2 + \mathbf{R}_4 \cdot q_2^*) + \mathbf{D}_2 v + \mathbf{D}_3 \dot{v}] = 0 \quad (B.13a) $$

$$ \dot{q}_2 - j\lambda_2^I q_2 - \epsilon\rho_2 q_2 = 0 \quad (B.13b) $$

$$ \epsilon \mathbf{L}_3^T [\mathbf{D}_1 (\mathbf{R}_1 \cdot q_1 + \mathbf{R}_2 \cdot q_1^* + \mathbf{R}_3 \cdot q_2 + \mathbf{R}_4 \cdot q_2^*) + \mathbf{D}_2 v + \mathbf{D}_3 \dot{v}] = 0 \quad (B.13b) $$

$$ \ddot{v} + \mu \dot{v} + \alpha v^3 + $$

$$ (\eta_1 + \eta_a)(R_{51}q_1 + R_{51}^*q_1^* + R_{53}q_2 + R_{53}^*q_2^*) - $$

$$ \eta_a(R_{71}q_1 + R_{71}^*q_1^* + R_{73}q_2 + R_{73}^*q_2^*) + $$

$$ (1 + \Omega_a^2)(R_{11}q_1 + R_{11}^*q_1^* + R_{13}q_2 + R_{13}^*q_2^*) + $$

$$ \gamma\Omega_2^2(R_{21}q_1 + R_{21}^*q_1^* + R_{23}q_2 + R_{23}^*q_2^*) - $$

$$ \Omega_a^2(R_{31}q_1 + R_{31}^*q_1^* + R_{33}q_2 + R_{33}^*q_2^*) + $$

$$ \epsilon [(\eta_1 + \eta_a + \mu)\dot{v} + (1 + \Omega_a^2)v + \alpha v^3 - $$

$$ (\eta_1 + \eta_a)(R_{51}q_1 + R_{51}^*q_1^* + R_{53}q_2 + R_{53}^*q_2^*) - $$

$$ (1 + \Omega_a^2)(R_{11}q_1 + R_{11}^*q_1^* + R_{13}q_2 + R_{13}^*q_2^*)] = 0 \quad (B.13c) $$

where $R_{ij}$ denotes the $(i, j)$ component of the matrix $\mathbf{R}$ and the vector $\mathbf{R}_i$. (respectively $\mathbf{L}_i$) denotes the $i^{th}$ column of the matrix $\mathbf{R}^T$ (respectively $\mathbf{L}^T$). Eq. (B.13) involves two complex variables $q_1$ and $q_2$ and one real variable $v$.

## B.3 ASYMPTOTIC ANALYSIS

The objective of this section is to characterize analytically the response regimes of the reduced model (B.13) in the simultaneous vicinity of the unstable natural frequencies $\lambda_1^I$ and $\lambda_2^I$. The proposed method which uses $\epsilon$ as a small parameter includes two steps. The first step (Sect. B.3.1) combines complexification and averaging methods to reduce the dynamics of the nonlinear system (B.13) by removing the oscillations with respect to the frequencies $\lambda_1^I$ and $\lambda_2^I$. The resulting dynamical system appears as a slow-fast system where the fast dynamics corresponds to the current time $t$ whereas the slow dynamics is related to the time $\epsilon t$. Note that the dynamics with respect to the frequencies $\lambda_1^I$ and $\lambda_2^I$ defines a third time scale denoted as the super fast dynamics. The reduced system can be analyzed using the geometric singular perturbation theory which is the objective of the second step (Sects. C.4.1 to B.3.5). That step begins with the definition of the critical manifold (Sect. C.4.1), which defines the space where the slow flow takes place. The associated stability properties are then considered in Sect. B.3.3. It continues with the characterization of the periodic or quasi-periodic regimes of the mechanical system (B.5) as stable fixed points of the slow flow (B.19) (Sect. B.3.4). Finally, in Sect. B.3.5, it is shown that the specific shape of the critical manifold enables the existence of more complicate steady-state responses named "relaxations oscillations", which exhibit strong modulated amplitudes on both frequency components.

### B.3.1 The complexification-averaging method and the slow-fast system

Our objective is to investigate the solution in the vicinity of two simultaneous natural frequencies $\lambda_1^I$ and $\lambda_2^I$ named $1: 1-1: 1$ resonances. We assume that $\lambda_1^I$ and $\lambda_2^I$ are incommensurable.

The complexification-averaging method [135] is used to make the slow dynamics of the system emerge from fast oscillating responses.

The complex modal variable $q_1$ (respectively $q_2$) may be written as a fast component with respect to the natural frequency $\lambda_1^I$ (respectively $\lambda_2^I$), and modulated by a slowly varying amplitude $\phi_1$ (respectively $\phi_2$),
that is

$$q_1 = \phi_1 e^{j\lambda_1^I t} \quad \text{and} \quad \dot{q}_1 = \dot{\phi}_1 e^{j\lambda_1^I t} + j\phi_1\lambda_1^I e^{j\lambda_1^I t} \qquad (B.14a)$$

$$q_2 = \phi_2 e^{j\lambda_2^I t} \quad \text{and} \quad \dot{q}_2 = \dot{\phi}_2 e^{j\lambda_2^I t} + j\phi_2\lambda_2^I e^{j\lambda_2^I t}. \qquad (B.14b)$$

As for the real variable $v$, new variables $v_1$ and $v_2$ are introduced as

$$v = v_1 + v_2 \qquad (B.15)$$

to capture frequency components with respect to $\lambda_1^I$ and $\lambda_2^I$ (fast components) respectively. The slowly changing complex amplitudes $\xi_1$ and $\xi_2$ are then respectively defined as

$$\xi_1 e^{j\lambda_1^I t} = \dot{v}_1 + j\lambda_1^I v_1 \qquad (B.16a)$$

$$\xi_2 e^{j\lambda_2^I t} = \dot{v}_2 + j\lambda_2^I v_2. \qquad (B.16b)$$

Inserting Eqs. (B.14), (B.15) and (B.16) into Eqs. (B.13) and averaging with respect to the frequencies $\lambda_1^I$ and $\lambda_2^I$ separately yields the following system of complex slow modulation with respect to $(\phi_1, \phi_2, \xi_1, \xi_2)$:

$$\dot{\phi}_1 = \epsilon F_1(\phi_1, \xi_1) \qquad (B.17a)$$

$$\dot{\phi}_2 = \epsilon F_2(\phi_2, \xi_2) \qquad (B.17b)$$

$$\dot{\xi}_1 = G_1(\phi_1, \xi_1, \xi_2, \epsilon) \qquad (B.17c)$$

$$\dot{\xi}_2 = G_2(\phi_2, \xi_2, \xi_1, \epsilon) \qquad (B.17d)$$

where the complex functions $F_1$, $F_2$, $G_1$ and $G_2$ are not explicitly given. Eqs. (B.17) characterize the dynamics of the slow complex amplitudes $\phi_1$, $\phi_2$, $\xi_1$ and $\xi_2$ of the physical variables $q_1$, $q_2$, $v_1$ and $v_2$ respectively.

By expanding Eqs. (B.17) using the following polar coordinates

$$\phi_1 = s_1 e^{j\delta_1} \qquad (B.18a)$$

$$\phi_2 = s_2 e^{j\delta_2} \qquad (B.18b)$$

$$\xi_1 = r_1 e^{j\theta_1} \qquad (B.18c)$$

$$\xi_2 = r_2 e^{j\theta_2}, \qquad (B.18d)$$

new equations of motion in terms of amplitudes $s_1, s_2, r_1$ and $r_2$ and phase differences $\Delta_1 = \delta_1 - \theta_1$ and $\Delta_2 = \delta_2 - \theta_2$ can be obtained as

$$\dot{s} = \epsilon \mathbf{f}(s, t) \qquad (B.19a)$$

$$\dot{t} = g(s, t, \epsilon) = g_0(s, t) + \epsilon g_1(s, t) \qquad (B.19b)$$

where $\mathbf{s} = (s_1, s_2)^T$ and $\mathbf{t} = (r_1, r_2, \Delta_1, \Delta_2)^T$. The functions $\mathbf{f} = (f_1, f_2)^T$, $\mathbf{g}_0 = (g_{01}, g_{02}, g_{03}, g_{04})^T$ and $\mathbf{g}_1 = (g_{11}, g_{21}, g_{31}, g_{41})^T$ take the following forms

$$f_1(s_1, r_1, \Delta_1) = a_1 s_1 + r_1(a_2 \cos \Delta_1 + a_3 \sin \Delta_1) \qquad (B.20)$$

$$f_2(s_2, r_2, \Delta_2) = a_4 s_2 + r_2(a_5 \cos \Delta_2 + a_6 \sin \Delta_2) \qquad (B.21)$$

$$g_{01}(s_1, r_1, \Delta_1) = H_{01}(r_1) + s_1(a_7 \cos \Delta_1 + a_8 \sin \Delta_1) \qquad (B.22)$$

$$g_{02}(s_2, r_2, \Delta_2) = H_{02}(r_2) + s_2(a_9 \cos \Delta_2 + a_{10} \sin \Delta_2) \qquad (B.23)$$

$$g_{03}(s_1, r_1, \Delta_1, r_2) = \frac{1}{r_1}(H_{03}(r_1, r_2) + s_1(a_8 \cos \Delta_1 - a_7 \sin \Delta_1)) \qquad (B.24)$$

$$g_{04}(s_2, r_2, \Delta_2, r_1) = \frac{1}{r_2}(H_{04}(r_1, r_2) + s_2(a_{10} \cos \Delta_2 - a_9 \sin \Delta_2)), \qquad (B.25)$$

with

$$H_{01}(r_1) = -\frac{\mu}{2}r_1, \qquad (B.26)$$

$$H_{02}(r_2) = -\frac{\mu}{2}r_2, \qquad (B.27)$$

$$H_{03}(r_1, r_2) = r_1 \left( \frac{\lambda_1^I}{2} - \frac{3\alpha}{8\lambda_1^{I3}} r_1^2 - \frac{3\alpha}{4\lambda_1^I \lambda_2^I} r_2^2 \right), \qquad (B.28)$$
$$H_{04}(r_1, r_2) = r_2 \left( \frac{\lambda_2^I}{2} - \frac{3\alpha}{4\lambda_1^{I2}\lambda_2^I} r_1^2 - \frac{3\alpha}{8\lambda_2^{I3}} r_2^2 \right). \qquad (B.29)$$

The coefficients $a_i$ for $i = 1, \dots, 10$, are reported in B.6 and the other functions are not given explicitly because they are not used in the following.

In the framework of the Geometric Singular Perturbation Theory, Eqs. (B.19) define a (4, 2)-fast-slow system where $t = (r_1, r_2, \Delta_1, \Delta_2)^T$ denote the fast variables and $s = (s_1, s_2)^T$ the slow ones. The slow variables are related to the modal components $q_1$ and $q_2$ whereas the fast variables are associated with the frequency components $v_1$ and $v_2$ of the NES variable $v$.

The system (B.19) corresponds to the fast time scale whereas for the slow time scale $\tau = \epsilon t$, the system reads

$$\mathbf{s}' = \mathbf{f}(\mathbf{s}, t) \qquad (B.30a)$$

$$\epsilon \mathbf{t}' = \mathbf{g}(\mathbf{s}, t, \epsilon) = \mathbf{g}_0(\mathbf{s}, t) + \epsilon \mathbf{g}_1(\mathbf{s}, t) \qquad (B.30b)$$

where $(\cdot)'$ denotes the derivative with respect to $\tau$. The slow and the fast subsystems are defined by considering $\epsilon = 0$, which yields for the slow subsystem

$$\mathbf{s}' = \mathbf{f}(\mathbf{s}, t) \qquad (B.31a)$$

$$\mathbf{0} = \mathbf{g}_0(\mathbf{s}, t), \qquad (B.31b)$$

which is a differential-algebraic equation, and for the fast subsystem

$$\dot{\mathbf{s}} = \mathbf{0} \qquad (B.32a)$$

$$\dot{\mathbf{t}} = \mathbf{g}_0(\mathbf{s}, t). \qquad (B.32b)$$

The slow and the fast subsystems will be complementary in the following analysis. Periodic or quasi-periodic steady-state regimes of the mechanical system (B.13) correspond to fixed points of the slow flow which are approximated by those of the slow subsystem. Indeed, if a fixed point of the slow flow with non-zero amplitude coordinates (i.e. $s_1, s_2, r_1$ and $r_2$) is reached, the resulting steady-state regime of the system (B.13) is approximated by a quasi-periodic regime containing the two harmonic components at the frequencies $\lambda_1^I$ and $\lambda_2^I$ respectively. If a fixed point with $s_1 = r_1 = 0$ (respectively $s_2 = r_2 = 0$) is reached, the resulting steady-state regime is approximated by a periodic regime containing only one harmonic component at the frequency $\lambda_2^I$ (respectively $\lambda_1^I$).

More complicate (amplitude and phase modulated) steady-state regimes of the mechanical system (B.13) will be explained, inter alia, by the existence of relaxation oscillations of the slow flow (see Sect. B.3.5), which are described combining slow and fast periods predicted by (E.30) and (E.31) respectively. In the case of only one unstable mode, those regimes are referred to as Strongly Modulated Responses (SMR) [188]. This term is kept in this work but it covers more situations which are described in Sect. B.3.5. Although SMR are also quasi-periodic regimes, we have chosen to reserve the "quasi-periodic" term for the situations evoked above in which a fixed point of the slow flow is reached.

### B.3.2 The critical manifold

The critical manifold is given as

$$C_0 := \{ (\mathbf{s}, \mathbf{t}) \in \mathbb{R}^2 \times \mathbb{R}^4 \mid \mathbf{g}_0(\mathbf{s}, \mathbf{t}) = 0 \}. \qquad (B.33)$$

It defines the domain where the slow flow takes place (see (E.30c)) and consists in the equilibrium of the fast flow (see (B.32b)). It plays a crucial role to characterize the response regimes under asymptotic assumptions.

The equations defining the critical manifold (E.32) can be simplified solving $g_{01} = 0$ and $g_{03} = 0$ (respectively $g_{02} = 0$ and $g_{04} = 0$) with respect to $\cos \Delta_1$ and $\sin \Delta_1$ (respectively $\cos \Delta_2$ and $\sin \Delta_2$), which gives

$$\cos \Delta_1 = - \frac{a_8 H_{03}(r_1, r_2) + a_7 H_{01}(r_1)}{s_1 (a_7^2 + a_8^2)} \quad \text{and} \quad \sin \Delta_1 = - \frac{a_8 H_{01}(r_1) - a_7 H_{03}(r_1, r_2)}{s_1 (a_7^2 + a_8^2)} \qquad (B.34)$$

(respectively

$$\cos \Delta_2 = - \frac{a_{10} H_{04}(r_1, r_2) + a_9 H_{02}(r_1)}{s_2 (a_9^2 + a_{10}^2)} \quad \text{and} \quad \sin \Delta_2 = - \frac{a_{10} H_{02}(r_1) - a_9 H_{04}(r_1, r_2)}{s_2 (a_9^2 + a_{10}^2)}. \qquad (B.35)$$
<figure><img src="image_74.png" /><figcaption>FIGURE B.2. The critical manifold in (a) the $(r_1, r_2, s_1)$-space (function (B.36a)) and (b) in the $(r_1, r_2, s_2)$-space (function (B.36b)). The set of parameters used are given in Sect. B.4.3, Eq. (B.58).</figcaption></figure>

Finally, compatibility equations relating $s_1$ and $s_2$ with $r_1$ and $r_2$ can be deduced from (B.34) and (B.35) as

$$s_1^2 = \frac{H_1(r_1, r_2)}{a_7^2 + a_8^2} \qquad (B.36a)$$

$$s_2^2 = \frac{H_2(r_1, r_2)}{a_9^2 + a_{10}^2} \qquad (B.36b)$$

where

$$H_1(r_1, r_2) = H_{01}^2(r_1) + H_{03}^2(r_1, r_2) \qquad (B.37a)$$

$$H_2(r_1, r_2) = H_{02}^2(r_2) + H_{04}^2(r_1, r_2). \qquad (B.37b)$$

The critical manifold $C_0$ relates the variables $s_1, s_2, \Delta_1$ and $\Delta_2$ to $r_1$ and $r_2$ (see Eqs. (B.34) to (B.37)) and appears as a two-dimensional parametric surface with respect to the two fast variables $r_1$ and $r_2$ which evolve in the six-dimensional variable space. It is hence possible to plot $C_0$ as two surfaces by projecting it on the three-dimensional spaces $(r_1, r_2, s_1)$ and $(r_1, r_2, s_2)$ using Eq. (B.37a) and Eq. (B.37b) respectively. Figure B.2 shows the corresponding two surfaces for given nominal parameter values of the considered system. The surface in the three-dimensional space $(r_1, r_2, s_1)$ (see Fig. B.2(a)) depicts $C_0$ with respect to the unstable mode $\lambda_1$. Note that, as observed when only one mode is unstable [a14], $C_0$ restricted to $r_2 = 0$ is S-shaped. Similar comments can be made considering the surface in the three-dimensional space $(r_1, r_2, s_2)$ related to the unstable mode $\lambda_2$. The presence of an S-shaped critical manifold will be used to show that SMR steady-state responses of the mechanical system can exist (as in the case of only one unstable mode, see for example [a14]).

Finally, the slow dynamics described by the slow subsystem can be only written in the $(r_1, r_2)$-plane, as performed hereafter in Sect. B.3.4.

### B.3.3 Stability of critical manifold

The points constituting $C_0$ are equilibrium point of the fast subsystem (E.31) and the stability properties can be analyzed from the Jacobian matrix of $g_0(s, t)$ with respect to $t$, denoted as $D_t g_0(s, t)$, and evaluated
at each point of $C_0$ as

$$ D_t g_0(\mathbf{s}, \mathbf{t}) = \begin{pmatrix} \frac{dH_{01}}{dr_1}(r_1) & -H_{03}(r_1, r_2) & 0 & 0 \\ 0 & 0 & \frac{dH_{02}}{dr_2}(r_2) & -H_{04}(r_1, r_2) \\ \frac{1}{r_1} \frac{\partial H_{03}}{\partial r_1}(r_1, r_2) & \frac{H_{01}(r_1)}{r_1} & \frac{1}{r_1} \frac{\partial H_{03}}{\partial r_2}(r_1, r_2) & 0 \\ \frac{1}{r_2} \frac{\partial H_{04}}{\partial r_1}(r_1, r_2) & 0 & \frac{1}{r_2} \frac{\partial H_{04}}{\partial r_2}(r_1, r_2) & \frac{H_{02}}{r_2}(r_2) \end{pmatrix} . \quad (B.38) $$

The critical manifold $C_0$ contains both regular points defined as $\det(D_t g_0(\mathbf{s}, \mathbf{t})) \neq 0$ and singular points characterized by $\det(D_t g_0(\mathbf{s}, \mathbf{t})) = 0$. Depending on the eigenvalues of $D_t g_0(\mathbf{s}, \mathbf{t})$, the set of regular points can be split into several parts (see definitions in Chap. 3 of [40]), including normally hyperbolic attracting parts (all the eigenvalues have strictly negative real parts), normally hyperbolic repelling parts (all the eigenvalues have strictly positive real parts), normally hyperbolic saddle type parts (neither attracting nor repelling) and non normally hyperbolic parts (at least one eigenvalue has a zero real part). Finally, the singular points for which the rank of $D_t g_0(\mathbf{s}, \mathbf{t})$ is equal to 3 (= 4 - 1), (the simplest possible rank deficiency) are denoted as fold points.

In general normally hyperbolic parts of $C_0$ constitute surfaces in the $(r_1, r_2)$-plane whereas points that are not normally hyperbolic (i.e. regular non normally hyperbolic points and singular points) are not isolated but form a one-dimensional curve inside $C_0$.

Fig. B.3 shows, for the same nominal parameter values of the considered system as above, the different parts of $C_0$ in the $(r_1, r_2)$-plane as

$$ C_0 = C_{0,r}^a \cup (F^1 \cup C_{0,r}^{st,1}) \cup (F^2 \cup C_{0,r}^{st,2}) \cup (F^3 \cup C_{0,r}^{st,3}) \cup (N_r \cup C_{0,r}^{st,4}) \quad (B.39) $$

where $C_{0,r}^a$ (Fig. B.3 - white background) is a (regular) normally hyperbolic attracting part and $C_{0,r}^{st,1}$ (respectively $C_{0,r}^{st,2}$ and $C_{0,r}^{st,3}$) (Fig. B.3 - salmon background) is a (regular) normally hyperbolic saddle type part that meets $C_{0,r}^a$ at the fold curve $F^1$ (respectively $F^2$ and $F^3$) (Fig. B.3 - dashed and continuous red lines). $C_{0,r}^{st,4}$ (Fig. B.3 - green background) is a normally hyperbolic saddle type part that meets $C_{0,r}^a$ at a (regular) non normally hyperbolic curve $N_r$ (Fig. B.3 - green line).

### B.3.4 Fixed points of the slow-flow in $C_0$

As mentioned previously, the fixed points of the slow flow (B.19) characterize periodic or quasi-periodic regimes of the mechanical system (B.5). Since $0 < \epsilon \ll 1$, these fixed points can be approximated by those of the slow subsystem (E.30) which are computed in this section.

The critical manifold $C_0$ relates the variables $s_1, s_2, \Delta_1$ and $\Delta_2$ to the variable $r_1$ and $r_2$. Hence, substituting Eq. (B.36) into (E.30), the dynamics of the slow subsystem can be considered only with respect to the variables $r_1$ and $r_2$ as

$$ D_r \mathbf{H}(\mathbf{r}) \mathbf{r}' = \mathcal{F}(\mathbf{r}) \quad (B.40) $$

where $\mathbf{r} = (r_1, r_2)^T$,

$$ D_r \mathbf{H}(\mathbf{r}) = \begin{pmatrix} \frac{\partial H_1}{\partial r_1}(r_1, r_2) & \frac{\partial H_1}{\partial r_2}(r_1, r_2) \\ \frac{\partial H_2}{\partial r_1}(r_1, r_2) & \frac{\partial H_2}{\partial r_2}(r_1, r_2) \end{pmatrix} \quad (B.41) $$

is the Jacobian matrix of the function vector $\mathbf{H}(\mathbf{r}) = (H_1(\mathbf{r}), H_2(\mathbf{r}))^T$ and

$$ \mathcal{F}(\mathbf{r}) = \begin{pmatrix} 2\sqrt{H_1(r_1, r_2)(a_7^2 + a_8^2)} \\ 2\sqrt{H_2(r_1, r_2)(a_9^2 + a_{10}^2)} \end{pmatrix} \mathbf{h}(\mathbf{r}). \quad (B.42) $$

In the above equation the vector function $\mathbf{h}(\mathbf{r})$ is obtained by replacing in $\mathbf{f}(\mathbf{s}, \mathbf{r})$ the variables $s_1, s_2, \Delta_1$ and $\Delta_2$ by their expressions as functions of $r_1$ and $r_2$.

From Eq. (B.40), it is possible to detect fixed points or equivalently regular fixed points as

$$ \mathcal{F}(\mathbf{r}) = 0, \quad (B.43a) $$

$$ \det(D_r \mathbf{H}(\mathbf{r})) \neq 0. \quad (B.43b) $$
<figure><img src="image_75.png" /><figcaption>**FIGURE B.3.** The different parts of the critical manifold $C_0$ are represented in the $(r_1, r_2)$-plane, as defined by (B.39) : $C_{0,r}^a$ (white background) is a (regular) normally hyperbolic attracting part; $C_{0,r}^{st,1}$ (respectively $C_{0,r}^{st,2}$ and $C_{0,r}^{st,3}$) (salmon background) is a (regular) normally hyperbolic saddle type part that meets $C_{0,r}^a$ at a fold curve $F^1$ (respectively $F^2$ and $F^3$) (dashed and continuous red lines); $C_{0,r}^{st,4}$ (green background) is also a normally hyperbolic saddle type part which meets $C_{0,r}^a$ at the (regular) non normally hyperbolic curve $N_r$ (green line). The blue dashed lines $A_1^1$ (respectively $A_1^2$ and $A_1^3$) are the arrival curves associated with the red dashed lines $F_1^1$ (respectively $F_1^2$ and $F_1^3$) and the blue solid lines $A_2^1$ (respectively $A_2^2$) are the arrival curves associated with the red solid lines $F_2^1$ (respectively $F_2^2$). The stream plot of the vector field $D_rH^{-1}\mathcal{F}(r)$ is depicted with gray arrows. Blue dots represent the regular fixed points of the slow flow as defined by (B.43) (here unstable). The red dots are the points where the transversality condition (B.48) does not hold and the orange dot is the point where $\text{adj}(D_rH)\mathcal{F}(r) = 0$. Same parameters as in Fig. B.2.</figcaption></figure>

For the given nominal parameter values, the stream plot of the vector field $D_rH^{-1}\mathcal{F}(r)$ associated with the dynamical system (B.40) is depicted in Fig. B.3 by gray arrows and the regular fixed points are reported as blue dots. A stability analysis has been performed given that all the fixed points are unstable. In the following magenta dots will be used to represent stable fixed points. A stable fixed point of the slow flow (B.19) characterizes a periodic regime of the initial system (B.5) if one of the variables $r_i$ ($i = 1, 2$) is equal to zero; otherwise the regime is quasi-periodic with two frequency components.

Fixed points no longer exist in situations where

$$ \mathcal{F}(r) = 0, \qquad (B.44a) $$

$$ \det(D_rH(r)) = 0. \qquad (B.44b) $$

The points in which Eqs. (B.44) hold are called *folded singularities* in the literature of slow-fast dynamical systems. Considering the relation

$$ \det(D_rH(r)) = -4r_1r_2 \det(D_tg_0(s, t)), \qquad (B.45) $$

it follows that fold curves can be defined from $\det(D_rH(r)) = 0$ and that folded singularities correspond to situations in which fixed points and fold points coincide. Folded singularities are hints of particular solutions of such slow-fast systems, called *canard cycles* [16]. These solutions are not investigated in this paper.

### B.3.5 Relaxation oscillations

As suggested by the particular shape of the critical manifold, we examine in this section if specific solutions of the slow-fast system named relaxation oscillations can occur. These solutions characterize complicated motions named SMR of the mechanical system (B.5). We focus here on the limit case $\epsilon \to 0$. A relaxation oscillation is a solution consisting of successive fast and slow segments forming a closed loop in the $(r_1, r_2)$-plane
representation of the critical manifold. The fast segments satisfy the fast subsystem (E.31) whereas the slow segments are obtained from the slow subsystem (E.30). Hence the slow segments take place on $C_0$ and the fast ones are outside of $C_0$. A switch from a slow segment to a fast one can occur when the trajectory crosses a fold line. In this case the fast segment appears as a jump from a fold point (in $C_0$) to an (arrival) point in $C_0$.

### B.3.5.1 Transversality condition

The transversality condition is defined as a crossing condition with respect to a fold curve $F^i$ ($i = 1, 2$ or $3$) of the slow dynamics. It is obtained writing that the slow flow, projected in the $(r_1, r_2)$-plane, is not tangent to a fold curve $F^i$ ($i = 1, 2$ or $3$) (i.e. a jump can occur).

The slow subsystem (B.40) being by definition singular on the fold curve, it is first desingularized introducing the time rescaling $\tau_1 = \tau/\det(D_rH(r))$, which leads to

$$ \mathbf{r}' = \operatorname{adj}(D_r(\mathbf{r})\mathbf{H}) \mathcal{F}(\mathbf{r}) \qquad (B.46) $$

where $\operatorname{adj}(.)$ denotes the adjoint matrix of $(.)$ and $(.)'$ is now the derivative with respect to the new time $\tau_1$.

Recalling that a fold curve in the $(r_1, r_2)$-plane is defined as $\det(D_rH(r)) = 0$ (see Sect. B.3.4), the transversality condition at a point $\mathbf{r}$ on a fold curve $F^i$ ($i = 1, 2$ or $3$) is written as follows

$$ \operatorname{grad}_{\mathbf{r}}(\det(D_{\mathbf{r}}\mathbf{H}(\mathbf{r}))) \cdot [\operatorname{adj}(D_{\mathbf{r}}\mathbf{H}) \mathcal{F}(\mathbf{r})] \neq 0 \qquad (B.47) $$

where $\operatorname{grad}_{\mathbf{r}}(.)$ denotes the gradient of $(.)$ with respect to $\mathbf{r}$ and $\cdot$ denotes the usual dot product in $\mathbb{R}^2$. The points where the transversality condition does not hold are obtained solving

$$ \operatorname{grad}_{\mathbf{r}}(\det(D_{\mathbf{r}}\mathbf{H}(\mathbf{r}))) \cdot [\operatorname{adj}(D_{\mathbf{r}}\mathbf{H}) \mathcal{F}(\mathbf{r})] = 0 \qquad (B.48) $$

for $\operatorname{adj}(D_rH) \mathcal{F}(\mathbf{r}) \neq 0$. The points defined by (B.48) are depicted in red on Fig. B.3 for the given nominal parameter values. These points are isolated and divide the fold curves into two parts

$$ F^i = F_1^i \cup F_2^i. \qquad (B.49) $$

In Fig. B.3 the red dashed (respectively solid) lines are arbitrarily chosen to depict $F_1^i$, $F_2^i$ and $F_3^i$ (respectively $F_1^D$, $F_2^D$ and $F_3^D$). Moreover, we can see that the stream plot of the vector field $D_rH^{-1}\mathcal{F}(\mathbf{r})$ can head towards or against the fold lines.

If $\operatorname{adj}(D_rH) \mathcal{F}(\mathbf{r}) = 0$, the considered points are no longer points where the transversality condition does not hold, but fixed points of the desingularized system (B.46), which are also hints of possible *canard cycles*. In Fig. B.3, a point, depicted as an orange dot, corresponds to this situation.

### B.3.5.2 Jumps and arrival curves

When the slow dynamics described by (B.40) reaches a point $\mathbf{r}^J = (r_1^J, r_2^J)^T$ of a part $F_j^i$ (with $i = 1, 2, 3$ and $j = 1, 2$) of the fold curve $F^i$ (this can occur if the stream plot of the vector field $D_rH^{-1}\mathcal{F}(\mathbf{r})$ heads towards $F_j^i$), the trajectory can leave the critical manifold $C_0$ (and the slow epoch) and undergoes a jump (a fast episode) before reaching another part of $C_{0,r}^a$. The point $\mathbf{r}^J$ is named *jump point*. The jump is described by the fast subsystem (E.31) with the jump point $\mathbf{r}^J$ as initial condition. Hence the possible drop point, denoted as $\mathbf{r}^D = (r_1^D, r_2^D)^T$, of the fast dynamics on $C_{0,r}^a$ is determined as a fixed point of the fast subsystem (E.31), that is

$$ \mathbf{g}_0(\mathbf{s}^J, \mathbf{t}) = 0 \qquad (B.50) $$

where $\mathbf{s}^J$ is obtained from Eqs. (B.36) evaluated at $\mathbf{r}^J$.

Manipulating Eq. (B.50), it can be shown that the drop point $\mathbf{r}^D$ is defined as the root of the following set of equations

$$ H_1(r_1^D, r_2^D) = H_1(r_1^J, r_2^J) \qquad (B.51a) $$

$$ H_2(r_1^D, r_2^D) = H_2(r_1^J, r_2^J), \qquad (B.51b) $$

where the functions $H_1$ and $H_2$ are defined by Eqs. (B.37). The solutions obtained solving Eqs. (B.51) for each point of $F_j^i$ define the *arrival curve* denoted as $A_j^i$. The arrival curves are also reported in Fig. B.3. The blue dashed line $A_1^1$ (respectively $A_1^2$ and $A_1^3$) is the arrival curve associated with the red dashed line $F_1^1$ (respectively $F_1^2$ and $F_1^3$) and the blue solid line $A_2^1$ (respectively $A_2^2$) is the arrival curve associated with the red solid line $F_2^1$ (respectively $F_2^2$).

The arrival curve (or a part of it) is stable, and therefore potentially reached by the trajectory at the end of a fast jump, if it is located on $C_{0,r}^a$, and unstable if it is located on $C_{0,r}^{st,i}$ ($i = 1, 2, 3$ or $4$).
### B.3.5.3 Description of a possible relaxation oscillation scenario

The previous results, leading to the general form of the critical manifold depicted in Figs. B.2 and B.3, allow us to predict different possible scenarios of relaxation oscillations. The latter can be defined from elementary relaxation oscillations.

An elementary relaxation oscillation is defined as a succession of fast, slow, fast and slow epochs starting from and dropping on the same part $F_j^i$. This can for example occur in Fig. B.3, considering the case where the slow dynamics in $C_{0,r}^a$ leads the trajectory to a jump point $P_1$ of $F_1^1$. Next a fast jump occurs before reaching a drop point $P_2$ of $A_1^1$. Then, from $P_2$, the slow dynamics reaches the fold curve $F_2^1$ at the jump point $P_3$, followed by a jump to a drop point $P_4$ of $A_2^1$. Finally the slow dynamics reaches the fold curve $F_1^1$ at a jump point $P_5$. If $P_1 = P_5$, the elementary relaxation oscillation defines a (sustained) simple relaxation oscillation.

In the case of Fig. B.3, it is easy to verify that elementary relaxation oscillations can occur in the neighborhood of $C_{0,r}^{st,1}$ and $C_{0,r}^{st,2}$, but not in the neighborhood of $C_{0,r}^{st,3}$. Indeed, a necessary condition to observe relaxation oscillations on both sides of a normally hyperbolic saddle type part $C_{0,r}^{st,i}$ ($i = 1, 2, 3$) of the critical manifold is that both $F_1^i$ and $F_2^i$ contain parts where the vector field $D_rH^{-1}\mathcal{F}(r)$ is incoming.

Unlike the case where the slow dynamics can be described in a one-dimensional space, here the critical manifold $C_0$ appears as a two-dimensional parametric surface. Consequently the condition $P_1 = P_5$ is not always satisfied. Moreover a simple relaxation oscillation can appear after a transient phase including successive elementary relaxation oscillations. We have not developed a method to characterize the existence of such simple relaxation oscillations.

More complex (sustained) relaxation oscillation scenarios may be considered, such as a double relaxation oscillation which is defined from two successive elementary relaxation oscillations occurring alternatively in the neighborhoods of two different $C_{0,r}^{st,i}$.

The next section provides an overview of the possible steady-states regimes by means of numerical simulations and aims at understanding them with regard to the theoretical results presented in this section. Note that basins of attraction of the possible steady-states regimes must be considered besides the previous local analysis. Indeed, for a given set of parameters, the system may have several stable attractors (stable fixed points or sustained relaxation oscillations) and the initial conditions determine which one is reached.

## B.4 NUMERICAL RESULTS

In this part, numerical simulations are performed both to validate and to complete the asymptotic analysis presented in Sect. B.3 by means of direct numerical integration of the slow flow (B.19) and of the initial system (B.12). For the numerical simulations, we fix $\epsilon = 0.01$ and the other parameters are always chosen to ensure that the first two modes of the primary system (B.3a-B.3d) undergo Hopf bifurcations at $\gamma = \gamma_{Hopf,1}$ and $\gamma = \gamma_{Hopf,2}$ respectively (i.e. the real parts $\rho_1$ and $\rho_2$ (see Eq. (B.9)) of the eigenvalues become greater than zero at $\gamma_{Hopf,1}$ and $\gamma_{Hopf,2}$ respectively). The numerical values of $\gamma_{Hopf,1}$ and $\gamma_{Hopf,2}$ are given for each situation under study in the following sections.

In all the examples presented below, the friction coefficient $\gamma$ is chosen such that $\gamma_{Hopf,1} < \gamma_{Hopf,2} < \gamma$. This choice validates the use of the reduced model (B.13) and the slow-flow formulation where $s_1$ and $s_2$ represent the two unstable modes of the primary system and $r_1$ and $r_2$ represent the frequency components (related to the two unstable modes of the primary system) of the relative displacement between the NES and the primary system.

Note that the primary structure (without NES) is linear, consequently it undergoes, when it is unstable, an exponential growth (i.e. the solution is not bounded). Therefore, all the responses of the coupled system shown in this section, which show bounded amplitudes, correspond to a quench of the "double-instability" by the single NES. These bounded motions can occur through periodic responses, quasi-periodic responses or complex SMR. Another type of response is shown at the end of the section (Sect B.4.4).

Sect. B.4.1 investigates the competition between two stable fixed points of the slow flow corresponding to the competition between two periodic or quasi-periodic regimes for the mechanical system; Sect. B.4.2 analyzes the competition between a sustained relaxation oscillation and a stable fixed point of the slow flow corresponding to the competition between an SMR and a periodic or quasi-periodic regime; in Sect. B.4.3 the competition between two kinds of sustained relaxation oscillations, corresponding to the competition between two SMR, is addressed. Finally, the case of the bifurcation of a fixed point to a limit cycle of the slow subsystem, corresponding to another type of amplitude and phase modulated response of the mechanical system, is presented in Sect. B.4.4.
### B.4.1 Competition between two stable fixed points of the slow flow

In this section the following set of parameters is used:

<table><tr><td>η<sub>1</sub> = 0.01</td><td>η<sub>2</sub> = 0.07</td><td>η<sub>3</sub> = 0.01</td><td>η<sub>4</sub> = 0.07</td><td>(B.52a)</td></tr><tr><td>Ω<sub>2</sub> = 0.93</td><td>Ω<sub>3</sub> = 1</td><td>Ω<sub>4</sub> = 1.7</td><td></td><td>(B.52b)</td></tr><tr><td>η<sub>a</sub> = 0.01</td><td>η<sub>a</sub> = 0.95</td><td></td><td></td><td>(B.52c)</td></tr><tr><td>σ = 1</td><td>γ = 0.09075</td><td></td><td></td><td>(B.52d)</td></tr><tr><td>μ = 0.395</td><td>α = 7.</td><td></td><td></td><td>(B.52e)</td></tr></table>

For this set of parameters, $γ_{Hopf,1} = 0.081$, $γ_{Hopf,2} = 0.088$ and the invariant manifold $C_0$ admits a decomposition in the form of (B.39). The fixed points of the slow flow are obtained solving Eq. (B.43), which gives four unstable fixed points including the trivial point (0,0), and two stable fixed points, denoted as $\mathbf{r}_{s,1}^*$ and $\mathbf{r}_{s,2}^*$, whose coordinates in the $(r_1, r_2)$-plane are respectively

$$ \mathbf{r}_{s,1}^* = (0.26, 0) \quad \text{and} \quad \mathbf{r}_{s,2}^* = (0, 0.45). \qquad (B.53) $$

The two stable fixed points of the slow flow correspond to a mitigation of the two unstable modes of the primary system. Note that if $\mathbf{r}_{s,1}^*$ is reached, the first unstable mode is mitigated through a periodic response of the physical system (B.12)<fnref n="1" /> whereas the second unstable mode is completely suppressed. If $\mathbf{r}_{s,2}^*$ is reached, the first unstable mode is completely suppressed whereas the second unstable mode is mitigated through a periodic response. The associated critical manifold $C_0$ is depicted on Fig. B.4(b) and Fig. B.5(b). The stable fixed points are represented by magenta dots and the unstable fixed points by blue dots. Two stable and one unstable fixed points are in the neighborhood of the trivial point (0,0). These four fixed points are on the (regular) normally hyperbolic attracting part $C_{0,r}^a$. The other two fixed points are far from (0,0). The stream plot of the vector field of Eq. (B.40), also reported in Fig. B.4(b) and Fig. B.5(b), suggests that each stable fixed point admits a basin of attraction in the neighborhood of the trivial point (0,0).

The direct numerical integration of the slow flow (B.19) is performed to localize, on the critical manifold $C_0$, the basin of attraction of each stable fixed point.

The time series of $s_1$, $s_2$, $r_1$ and $r_2$ obtained from two different initial conditions $\mathbf{r}_{0,1} = (r_1(0), r_2(0)) = (0.02, 0.1)$ and $\mathbf{r}_{0,2} = (0.02, 0.2)$ are reported respectively on Fig. B.4(a) and Fig. B.5(a), and the corresponding trajectories in the $(r_1, r_2)$-plane are plotted on Fig. B.4(b) and Fig. B.5(b). In the latter, the initial condition points $\mathbf{r}_{0,1}$ and $\mathbf{r}_{0,2}$ are represented by light blue dots and the trajectories of the slow flow are depicted using a color gradient from light blue at $t=0$ to black at $t=30000$. These trajectories give a good idea of the basins of attraction and of the role played by the unstable fixed points. The basin of attraction of $\mathbf{r}_{s,1}^*$ is located between the $r_1$-axis and the streamline that links the trivial point to the unstable fixed point located at (0.16, 0.25); the basin of attraction of $\mathbf{r}_{s,2}^*$ is located between the same streamline and the $r_2$-axis.

Finally, the time series obtained solving the complete equations (B.12) are also reported in Fig. B.4(a) and Fig. B.5(a) following the remarks below.

From Eqs. (B.14) to (B.16) and (B.18), the variables of the complete model (B.12) can be expressed as functions of the slow-flow variables as follows

$$ q_1 = s_1 e^{j(\lambda_1^I t + \delta_1)}, \qquad q_2 = s_2 e^{j(\lambda_2^I t + \delta_2)}, \qquad (B.54) $$

$$ v_1 = \frac{r_1}{\lambda_1^I} \sin(\lambda_1^I t + \theta_1) = A_{v_1} \sin(\lambda_1^I t + \theta_1), \qquad v_2 = \frac{r_2}{\lambda_2^I} \sin(\lambda_2^I t + \theta_2) = A_{v_2} \sin(\lambda_2^I t + \theta_2). \qquad (B.55) $$

To analyze the results, the modulus of the variables $q_1$ and $q_2$ (orange lines) of the system (B.12) must be compared to the variables $s_1$ and $s_2$ (blue lines) of the slow flow (B.19) respectively. Moreover, the spectrograms of the signal $v(t)$ (resulting from the numerical integration of Eq. (B.12)) are then computed and the amplitudes $A_{v_1}$ and $A_{v_2}$ of the components $v_1(t)$ and $v_2(t)$ related to the frequencies $\lambda_1^I$ and $\lambda_2^I$ respectively are extracted. The variables $\lambda_1^I A_{v_1}$ and $\lambda_2^I A_{v_2}$ (black dashed lines) are to be compared to the variables $r_1$ and $r_2$ (green lines) of the slow flow (B.19).

In each case the behavior of $s_1$ (respectively $r_1$, $s_2$ and $r_2$) is found close to the behavior of $q_1$ (respectively $\lambda_1^I A_{v_1}$, $q_2$ and $\lambda_2^I A_{v_2}$), which allows us to validate the theoretical development presented in Section B.3.

Even if the previous comparisons show a good agreement, quantitative differences are nevertheless observed. We recall that two steps separate the initial mechanical system (B.12) and its approximated slow

<footnote n="1">1. Indeed, if $r_1 = 0$ or $r_2 = 0$, the non averaged initial system undergoes periodic oscillations.</footnote>
<figure><img src="image_76.png" /><figcaption>**FIGURE B.4.** Competition between two stable fixed points of the slow flow (1st example). (a) Amplitudes $r_1, s_1$ (top) and $r_2, s_2$ (bottom) as functions of time with the initial condition $\mathbf{r}_{0,1} = (0.02, 0.1)$ obtained by solving numerically the slow flow (B.19). The time series of $q_1$ (top) and $q_2$ (bottom) obtained solving numerically Eqs. (B.12) are also shown. (b) Similar representation as in Fig. B.3 with the addition of the trajectories (solid lines with a color gradient from light blue at $t = 0$ to black at $t = 30000$) shown in (a); the magenta dots are the stable fixed points of the slow flow. The parameters (G.45) are used.</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>**FIGURE B.5.** Competition between two stable fixed points of the slow flow (2nd example). Same caption as in Fig. B.4 with the initial condition $\mathbf{r}_{0,2} = (0.02, 0.2)$.</figcaption></figure>

flow (B.19) : the reduction of the dynamics keeping only the unstable modes and the averaging with respect to the frequencies of these two unstable modes. Additional numerical simulations (not displayed here) show that the quantitative differences observed here (and also in the following sections) are mainly due to the averaging step.

## B.4.2 Competition between a sustained relaxation oscillation and a stable fixed point of the slow flow

To illustrate the competition between sustained relaxation oscillations and a stable fixed point we use the following parameters

$$ \eta_1 = 0.01 \qquad \eta_2 = 0.07 \qquad \eta_3 = 0.01 \qquad \eta_4 = 0.07 \qquad (B.56a) $$
<table><tr><td>Ω<sub>2</sub> = 0.93</td><td>Ω<sub>3</sub> = 1</td><td>Ω<sub>4</sub> = 1.7</td><td>(B.56b)</td></tr><tr><td>η<sub>a</sub> = 0.01</td><td>Ω<sub>a</sub> = 0.95</td><td></td><td>(B.56c)</td></tr><tr><td>σ = 1</td><td>γ = 0.09</td><td></td><td>(B.56d)</td></tr><tr><td>μ = 0.25</td><td>α = 7.</td><td></td><td>(B.56e)</td></tr></table>

The parameters of the primary system are here unchanged compared to those of Sect. B.4.1 (only the bifurcation parameter $γ$ and the NES parameters $μ$ and $α$ are modified). Therefore, the values of the Hopf bifurcation points are still $γ_{Hopf,1} = 0.081$ and $γ_{Hopf,2} = 0.088$. Moreover, the invariant manifold $C_0$ admits again a decomposition in the form of (B.39); the slow flow now admits six fixed points, among which only one is stable. The latter, denoted as $\mathbf{r}_{s,3}^*$, has the following coordinates in the $(r_1, r_2)$-plane :

$$ \mathbf{r}_{s,3}^* = (0, 0.57). \qquad (B.57) $$

If $\mathbf{r}_{s,3}^*$ is reached, only the first unstable mode is completely suppressed whereas the second unstable mode is mitigated through a periodic response. The critical manifold $C_0$ is reported in Fig. B.6(b) and Fig. B.7(b) where the stable fixed point is again represented by a magenta dot and the unstable points by blue dots. The stable fixed point and two unstable fixed points (excluding the trivial point) are in the neighborhood of the trivial point (0,0), the other unstable points being on the contrary far from (0,0). The stream plot of the vector field of Eq. (B.40) suggests that, depending on the initial condition point, the system can reach the stable fixed point $\mathbf{r}_{s,3}^*$ or initiate another type of oscillations.

Direct numerical integration of the slow flow (B.19) is performed to localize, on the critical manifold $C_0$, the basin of attraction of the stable fixed point and to characterize the other type of oscillations. The time series obtained from the slow flow (B.19) considering the initial condition $\mathbf{r}_{0,3} = (0.01, 0.15)$ are plotted in Fig. B.6(a) and, using a color gradient, the corresponding trajectory in the $(r_1, r_2)$-plane is depicted in Fig. B.6(b). It can be seen that the stable fixed point is reached.

Results related to a second initial condition point $\mathbf{r}_{0,4} = (0.01, 0.03)$ are reported in Fig. B.7. Again, the time series are plotted in Fig. B.7(a) and the trajectory in the $(r_1, r_2)$-plane in Fig. B.7(b). To better understand the observed behavior, three-dimensional plots are also used (see Figs. B.7(c) and B.7(d)), where the trajectories are plotted in the $(r_1, r_2, s_1)$-space and the $(r_1, r_2, s_2)$-space respectively and superimposed to the surfaces, defined by Eqs. (B.36), characterizing the invariant manifold $C_0$. Moreover, using the same color code as in Fig. B.3, the fold and arrival curves as well as the points where the transversality condition is lost are depicted in Figs. B.7(c) and B.7(d).

Fig. B.7(b) shows that after transient elementary relaxation oscillations (described in Sect. B.3.5.3) with respect to the normally hyperbolic saddle type part $C_{0,r}^{st,1}$ of the critical manifold, a sustained simple relaxation oscillation (described in Sect. B.3.5.3) occurs along the straight line $r_2 = 0$ ($s_2$ and $r_2$ tend towards zero as visible in Fig. B.7(a)). Observing Figs. B.7(c) and B.7(d), we can see that this simple relaxation oscillation is due to an S-shape of the critical manifold in the $(r_1, r_2, s_1)$-space and in the plane $r_2 = 0$. This means that the first unstable mode here is completely suppressed while the second unstable mode is mitigated through a sustained simple relaxation oscillation which corresponds to an SMR for the initial physical system (B.12).

As in the previous case studied in Sect. B.4.1, the boundary between the basins of attractions of the two attractors appears to be the connection between two unstable fixed points, i.e. the trivial fixed point and the fixed point located at (0.22, 0.18).

Here again, the predictions obtained with the slow flow (B.19) are in good agreement with the simulations of the complete system (B.12) (see Figs. B.6(a) and B.7(a)).

### B.4.3 Competition between two kinds of sustained relaxation oscillations

In this section the following parameters are used

<table><tr><td>η<sub>1</sub> = 0.01</td><td>η<sub>2</sub> = 0.07</td><td>η<sub>3</sub> = 0.01</td><td>η<sub>4</sub> = 0.07</td><td>(B.58a)</td></tr><tr><td>Ω<sub>2</sub> = 0.93</td><td>Ω<sub>3</sub> = 1</td><td>Ω<sub>4</sub> = 1.7</td><td></td><td>(B.58b)</td></tr><tr><td>η<sub>a</sub> = 0.01</td><td>Ω<sub>a</sub> = 0.95</td><td></td><td></td><td>(B.58c)</td></tr><tr><td>σ = 0.85</td><td>γ = 0.0935</td><td></td><td></td><td>(B.58d)</td></tr><tr><td>μ = 0.28</td><td>α = 5</td><td></td><td></td><td>(B.58e)</td></tr></table>

for which the slow flow (B.19) has six fixed points, all unstable. The values of the Hopf bifurcation points are $γ_{Hopf,1} = 0.078$ and $γ_{Hopf,2} = 0.088$. The critical manifold $C_0$ is reported in Fig. B.8(b) and Fig. B.9(b). Note that there are no unstable fixed points in the neighborhood of the trivial point. Obviously in this case, because
<figure><img src="image_78.png" /><figcaption>**FIGURE B.6.** Competition between a sustained relaxation oscillation and a stable fixed point of the slow flow (1st example). Same caption as in Fig. B.4 with the initial condition $\mathbf{r}_{0,3} = (0.01, 0.15)$. The set of parameters (F.4) is used.</figcaption></figure>

the slow flow has no stable fixed point, the instability mitigation cannot be performed through quasi-periodic or periodic regimes.

In practice, depending on the chosen initial conditions, two kinds of sustained relaxation oscillations can be observed by means of numerical simulations of the slow flow (B.19).

Fig. B.8 illustrates the first scenario. The time series are obtained from the initial condition $\mathbf{r}_{0,5} = (0.05, 0.05)$ and are reported in Fig. B.8(a). This scenario is similar to that observed in Sect. B.4.2 and reported in Fig. B.7: after transient elementary relaxation oscillations, the second unstable mode is completely suppressed (see Fig. B.8(a)) whereas the first unstable mode is mitigated through a sustained simple relaxation oscillation taking place along the straight line $r_2 = 0$ (see the black curve in Fig. B.8(c)). This simple relaxation oscillation is due to the S-shape of the critical manifold in the $(r_1, r_2, s_1)$-space and in the plane $r_2 = 0$. It consists of four parts: two fast jumps (from the fold curve $F_1^1$ to the arrival curve $A_1^1$ and from the fold curve $F_2^1$ to the arrival curve $A_2^1$) and two slow epochs on $C_{0,r}^a$.

Fig. B.9 illustrates the second scenario. The time series are obtained from the initial condition $\mathbf{r}_{0,6} = (0.1, 0.1)$ and are reported in Fig. B.9(a). After a quick transient phase, the two unstable modes are mitigated through complex relaxation oscillations with respect to the normally hyperbolic saddle type parts $C_{0,r}^{st,1}$ and $C_{0,r}^{st,2}$ of the critical manifold successively. The complex relaxation oscillations here consist of three parts: an elementary relaxation oscillation due to the S-shape of the critical manifold in the $(r_1, r_2, s_2)$-space (see Fig. B.9(d)), an elementary relaxation oscillation due to the S-shape of the critical manifold in the $(r_1, r_2, s_1)$-space (see Fig. B.9(c)) and transient bursting oscillations (see Figs. B.9(b) and B.9(c)). This sustained complex relaxation oscillation regime is due to successive targeted energy transfers, first from the first unstable mode to the NES (principally when the jumps from $F_1^2$ to $A_1^2$ and from $F_2^2$ to $A_2^2$ occur) and then from the second mode to the NES (when the jumps from $F_1^1$ to $A_1^1$ occur). Note that those energy transfers do not occur simultaneously as shown in Fig. B.9(a).

From [40], bursting oscillations are defined as patterns of time series alternating between near steady-state and rapid oscillatory phases. Here, the steady-state phases are the slow evolution parts on $C_{0,r}^a$ of the relaxation oscillation cycles and the rapid oscillatory phases correspond to the fast jumps which are now oscillating. This phenomenon is related to the normally hyperbolic saddle type part of $C^0$, $C_{0,r}^{st,4}$, that meets $C_{0,r}^a$ at the (regular) non normally hyperbolic curve $N_r$ (see Fig. B.9(b)). Indeed, when a part of the arrival curve is in this area, it is made up of unstable fixed points of the fast subsystem (E.31). In this case, the attractor of the fast dynamics during a jump, at constant values of $s_1$ and $s_2$ (see Eq. (B.32a)), is no longer an equilibrium point but a limit cycle oscillation (LCO). Then, when the jump is finished, because $s_1$ and $s_2$ vary again the system leaves $C_{0,r}^{st,4}$ and the LCO undergoes a bifurcation to a new steady-state (slow) phase. Here the bursting oscillations are transient but depending on the parameter values, they can also been sustained. This scenario has not been studied in this paper.

The boundary between the basins of the two attractors does not appear as clearly as in the previous situations where it corresponded to the connection between two unstable points. Indeed, as noticed above,
<figure><img src="image_79.png" /><figcaption>FIGURE B.7. Competition between a sustained relaxation oscillation and a stable fixed point of the slow flow (2nd example). (a) and (b) Same captions as in Fig. B.4 with the initial condition $r_{0,4} = (0.01, 0.03)$. In (c) and (d) the trajectory of the slow flow obtained from the numerical simulation of Eq. (B.19) is plotted in the ($r_1, r_2, s_1$)-space and the ($r_1, r_2, s_2$)-space respectively and superimposed to the surfaces given by Eqs. (B.37a) and (B.37b). The set of parameters (F.4) is used.</figcaption></figure>

there are no unstable fixed points in the neighborhood of the trivial point. Here, the boundary depends on the first reached point on the fold curve $F_1$ and on the direction of the transient relaxation oscillations along $C_{0,r}^{st,1}$, which can move downward (first scenario) or upward (second scenario).

Finally, here also, the predictions obtained with the slow flow (B.19) are in good agreement with the complete system (B.12) (see Figs. B.8(a) and B.9(a)).

### B.4.4 Bifurcation of a fixed point to a limit cycle of the slow subsystem

For this last example, the following parameters are used with two different values of the friction coefficient

$$ \eta_1 = 0.01 $$

$$ \Omega_2 = 0.93 $$

$$ \eta_2 = 0.07 $$

$$ \Omega_3 = 1 $$

$$ \eta_3 = 0.01 $$

$$ \Omega_4 = 1.7 $$

$$ \eta_4 = 0.07 $$

(B.59a)

(B.59b)
<figure><img src="image_80.png" /><figcaption>**FIGURE B.8.** Competition between two kinds of sustained relaxation oscillations (1st example). Same caption as in Fig. B.7 with the initial condition $\mathbf{r}_{0,5} = (0.05, 0.05)$. The set of parameters (B.58) is used.</figcaption></figure>

$$ \eta_a = 0.01 \qquad \Omega_a = 0.95 \qquad \sigma = 0.98 \qquad (B.59c) $$

$$ \gamma_1 = 0.09053 \text{ (1st example)} \qquad \gamma_2 = 0.09075 \text{ (2nd example)} \qquad (B.59d) $$

$$ \mu = 0.45 \qquad \alpha = 7. \qquad (B.59e) $$

With the parameters (B.59), the primary system undergoes Hopf bifurcations at $\gamma_{\text{Hopf},1} = 0.08$ and $\gamma_{\text{Hopf},2} = 0.086$. The associated critical manifolds $C_0$ are reported in Fig. B.10(b) and Fig. B.10(d) for $\gamma = 0.09053$ and $\gamma_2 = 0.09075$ respectively.

For $\gamma_1 = 0.09053$, the slow flow (B.19) admits a single stable fixed point having the following coordinates in the $(r_1, r_2)$-plane :

$$ \mathbf{r}_{s,4}^* = (0.06, 0.61), \qquad (B.60) $$

and five unstable fixed points. Note that if $\mathbf{r}_{s,4}^*$ is reached, both unstable modes are mitigated through a periodic response of the physical system which therefore undergoes a quasi-periodic regime even if the amplitude of the component corresponding to the first unstable mode is much smaller than that of the second unstable mode. Numerical simulations of the slow flow (B.19) corresponding to this situation are shown in Figs. B.10(a) and
<figure><img src="image_81.png" /><figcaption>**FIGURE B.9.** Competition between two kinds of sustained relaxation oscillations (2^nd^ example). Same caption as in Fig. B.7 with the initial condition r<sub>0,6</sub> = (0.05, 0.15). The set of parameters (B.58) is used.</figcaption></figure>

B.10(b).

Then, when $\gamma$ increases, the stable fixed point $\mathbf{r}_{s,4}^*$ of the slow flow undergoes a Hopf bifurcation and becomes unstable. For example, for $\gamma = 0.09075$, the slow flow has only unstable fixed points, four of which are in the neighborhood of (0, 0).

Numerical simulations of the slow flow (B.19) corresponding to this second case have been performed for different values of the initial conditions, showing that only one stable attractor is present whose basin of attraction includes a neighborhood of the trivial point (0,0) (as it was previously the case for the stable fixed point). The time series obtained from the initial condition $r_{0,7} = (0.05, 0.05)$ are reported in Figs. B.10(c) and Fig. B.10(d). The figures show that after the bifurcation of the stable fixed point (B.60), the slow flow undergoes an LCO (see the black curve in Fig. B.10(d)) which takes place in the attractive part $C_{0,r}^a$ of the critical manifold, and therefore only at the slow-time scale (there are no fast jumps, unlike for relaxation oscillations).
<figure><img src="image_82.png" /><figcaption>**FIGURE B.10.** *Bifurcation of a fixed point to a limit cycle of the slow subsystem.* Same caption as in Fig. B.4 with the initial condition $\mathbf{r}_{0,7} = (0.05, 0.05)$. The parameters (B.59) are used, with (a) and (b) $\gamma = 0.09053$ and (c) and (d) $\gamma = 0.09075$.</figcaption></figure>

Therefore, the LCO can be characterized theoretically only from the two-dimensional slow subsystem (B.40). This situation also corresponds to a modulated response, in terms of amplitude and phase, of the mechanical system (B.13), but here due to a limit cycle of the slow flow instead of relaxation oscillations of the slow flow as in the SMR case.

Finally, once again, the predictions obtained with the slow flow (B.19) are in good agreement with that of the complete system (B.12) (see Figs. B.10(a) and B.10(c)), which allows us to validate the asymptotic analysis presented in Sect. B.3.

## B.5 CONCLUSION

The possibility of mitigating simultaneously two unstable modes of a linear MDOF primary system using a nonlinear energy sink (NES) has been investigated in this paper by means of an asymptotic analysis complemented by numerical simulations. This study focuses on a 4-DOF phenomenological friction system having two unstable modes and coupled to a single NES.

The asymptotic analysis begins with the well-known biorthogonal transformation to diagonalize the primary system for the purpose of reducing its dynamics only to the two unstable modes. A complexification-averaging process is applied at the vicinity of the two 1:1-1:1 resonances at the natural frequencies of the unstable modes, leading to the slow flow of the system which takes the form of a (4,2)-fast-slow system. The
slow variables characterize the two unstable modes of the primary system whereas the four fast variables describe the relative motion between the NES and the primary system (amplitudes and phases of the two frequency components relative to the two unstable modes). The slow flow is finally analyzed by means of the geometric singular perturbation theory. The critical manifold $C_0$ is obtained and, contrary to the case of a single unstable mode, it consists in a 2-dimensional parametric surface with respect to the two fast variables $r_1$ and $r_2$ which evolve in the whole 6-dimensional variable space, and which characterize the frequency components of the NES dynamics with respect to the two unstable modes. The critical manifold appears as the union of a normally hyperbolic attracting domain and four normally hyperbolic saddle type domains (with finite dimensions). Two of them meet the attracting part at fold curves and it can be shown by analyzing, inter alia, the transversality condition, that these fold curves can be split into successive incoming and outgoing sections, enabling jumps and complex motions of the slow flow, such as successive slow and fast epochs (i.e. relaxation oscillations). Some relaxation oscillation scenarios have been deduced from the asymptotic analysis, among which simple relaxation oscillations around one of the normally hyperbolic saddle type domains, or multi-relaxation oscillations around two of the normally hyperbolic saddle type domains. These relaxation oscillation scenarios and the stable fixed points of the slow flow define possible sustained mitigated regimes of the physical system such as strongly modulated responses (SMR), complex SMR, periodic or quasi-periodic regimes.

Direct numerical integrations of the system have been performed. First of all, the simulations validated that it was possible to mitigate simultaneously two unstable modes by means of a single NES. The two unstable modes can be completely eliminated, or mitigated through one of the sustained regimes (periodic, quasi-periodic, SMR). Secondly, the simulations allowed us to study the notion of basin of attraction. Indeed, because of the two-dimensional nature of $C_0$ and the fact that, for a given set of parameters, the system may have several stable attractors (stable fixed points or sustained relaxation oscillations) this notion has to be considered. Three kinds of competitions have been observed, respectively between two stable fixed points of the slow flow, between a sustained relaxation oscillation and a stable fixed point of the slow flow and between two kinds of sustained relaxation oscillations. In each situation considered, a good agreement has been observed between the theoretical results and the numerical simulations, which validates the proposed asymptotic analysis. Moreover, the numerical simulations highlighted that the two unstable modes could also be mitigated through a strongly modulated response defined from a limit cycle of the slow system at the slow time scale.

As a general conclusion, when it is coupled to a linear system which has two unstable modes, a single NES is able to produce bounded motions of various kinds, thoroughly studied in this paper. In real life, unstable mechanical systems are in general nonlinear and therefore they undergo limit cycle oscillations which can be dangerous if they are of significant amplitude. As a perspective, from an engineering point of view, the NES may be designed to produce effective bounded motions, i.e. with a lower amplitude than that of the system without the NES. The present study might constitute a tool to perform this design. The amplitude of the oscillations of the resulting bounded motions may also be compared with those provided by other mitigation methods.

## ACKNOWLEDGEMENTS

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

## B.6 APPENDIX A. COEFFICIENTS OF THE SLOW-FOW SYSTEM (B.19)

The coefficients $a_i$ for $i = 1, \dots, 10$, characterizing the functions $f$ and $g$ which define the slow-fow system (B.19), are given by

$$a_1 = \frac{\eta_a L_{71}^I R_{51}^I}{\sigma} - \eta_a L_{51}^I R_{51}^I + \frac{\Omega_a^2 L_{71}^I R_{11}^I}{\sigma} - \Omega_a^2 L_{51}^I R_{11}^I - \frac{R_{51}^R \eta_a L_{71}^R}{\sigma} + R_{51}^R \eta_a L_{51}^R - \frac{R_{11}^R \Omega_a^2 L_{71}^R}{\sigma} + R_{11}^R \Omega_a^2 L_{51}^R + \gamma L_{61}^I R_{11}^I - \eta_1 L_{51}^I R_{51}^I - L_{51}^I R_{11}^I - \gamma R_{11}^R L_{61}^R + \eta_1 R_{51}^R L_{51}^R + R_{11}^R L_{51}^R + \rho_1, \quad (B.61)$$

$$a_2 = \frac{\Omega_a^2 L_{71}^I}{2\sigma\lambda_1^I} - \frac{\Omega_a^2 L_{51}^I}{2\lambda_1^I} + \frac{\eta_a L_{71}^R}{2\sigma} - \frac{1}{2}\eta_a L_{51}^R + \frac{\gamma L_{61}^I}{2\lambda_1^{-I}} - \frac{L_{51}^I}{2\lambda_1^I} - \frac{1}{2}\eta_1 L_{51}^R, \quad (B.62)$$

$$a_3 = \frac{\eta_a L_{71}^I}{2\sigma} - \frac{1}{2}\eta_a L_{51}^I - \frac{\Omega_a^2 L_{71}^R}{2\sigma\lambda_1^I} + \frac{\Omega_a^2 L_{51}^R}{2\lambda_1^I} - \frac{1}{2}\eta_1 L_{51}^I - \frac{\gamma L_{61}^R}{2\lambda_1^I} + \frac{L_{51}^R}{2\lambda_1^I}, \quad (B.63)$$
$$a_4 = \frac{\eta_a L_{73}^I R_{53}^I}{\sigma} - \eta_a L_{53}^I R_{53}^I + \frac{\Omega_a^2 L_{73}^I R_{13}^I}{\sigma} - \Omega_a^2 L_{53}^I R_{13}^I - \frac{R_{53}^R \eta_a L_{73}^R}{\sigma} + R_{53}^R \eta_a L_{53}^R - \frac{R_{13}^R \Omega_a^2 L_{73}^R}{\sigma} + R_{13}^R \Omega_a^2 L_{53}^R + \gamma L_{63}^I R_{13}^I - \eta_1 L_{53}^I R_{53}^I - L_{53}^I R_{13}^I - \gamma R_{13}^R L_{63}^R + \eta_1 R_{53}^R L_{53}^R + R_{13}^R L_{53}^R + \rho_2, \quad (B.64)$$

$$a_5 = \frac{\Omega_a^2 L_{73}^I}{2\sigma\lambda_2^I} - \frac{\Omega_a^2 L_{53}^I}{2\lambda_2^I} + \frac{\eta_a L_{73}^R}{2\sigma} - \frac{1}{2}\eta_a L_{53}^R + \frac{\gamma L_{63}^I}{2\lambda_2^I} - \frac{L_{53}^I - \frac{1}{2}\eta_1 L_{53}^R}{2\lambda_2^I}, \quad (B.65)$$

$$a_6 = \frac{\eta_a L_{73}^I}{2\sigma} - \frac{1}{2}\eta_a L_{53}^I - \frac{\Omega_a^2 L_{73}^R}{2\sigma\lambda_2^I} + \frac{\Omega_a^2 L_{53}^R}{2\lambda_2^I} - \frac{1}{2}\eta_1 L_{53}^I - \frac{\gamma L_{63}^R}{2\lambda_2^I} + \frac{\lambda_2^I L_{53}^R L_{53}^R}{2\lambda_2^I L_{53}^R}, \quad (B.66)$$

$$a_7 = -R_{51}^R \eta_a + R_{71}^R \eta_a - R_{11}^R \Omega_a^2 + R_{31}^R \Omega_a^2 - \gamma R_{21}^R \Omega_2^2 - \eta_1 R_{51}^R - R_{11}^R \quad (B.67)$$

$$a_8 = \eta_a R_{51}^I - \eta_a R_{71}^I + \Omega_a^2 R_{11}^I - \Omega_a^2 R_{31}^I + \gamma \Omega_2^2 R_{21}^I + \eta_1 R_{51}^I + R_{11}^I, \quad (B.68)$$

$$a_9 = -R_{53}^R \eta_a + R_{73}^R \eta_a - R_{13}^R \Omega_a^2 + R_{33}^R \Omega_a^2 - \gamma R_{23}^R \Omega_2^2 - \eta_1 R_{53}^R - R_{13}^R$$

and

$$a_{10} = \eta_a R_{53}^I - \eta_a R_{73}^I + \Omega_a^2 R_{13}^I - \Omega_a^2 R_{33}^I + \gamma \Omega_2^2 R_{23}^I + \eta_1 R_{53}^I + R_{13}^I. \quad (B.69)$$
*Submitted to Physica D : Nonlinear Phenomena (2024) [a13]*

# FAST-SLOW ANALYSIS OF PASSIVE MITIGATION OF SELF-SUSTAINED OSCILLATIONS BY MEANS OF A BISTABLE NONLINEAR ENERGY SINK

Baptiste BERGEOT and Sébastien BERGER

## ABSTRACT

This paper investigates the dynamic behavior of a Van der Pol oscillator (used as an archetypal self-sustained oscillator) coupled to a bistable nonlinear energy sink (BNES). We first show using numerical simulations that this system can undergo a multitude of motions including different types of periodic regimes and so-called strongly modulated responses (SMR) as well as chaotic regimes. We also show that a BNES can be much more efficient than a classical cubic NES but this is not robust since a little perturbation can switch the system from harmless to harmful situations. However, even in the most unfavorable cases, it is possible to find a set of parameters for which the BNES performs better than the NES.

A multiple time scales approach is then addressed to analyze the system. In this context, we show that the so-called Multiple Scale/Harmonic Balance Method (MSHBM) must be modified (compared to its usual use) to consider the specific feature of the BNES, i.e., that it can have a nonzero-mean oscillating motion. This allows us to derive a so-called amplitude-phase modulation dynamics (APMD) which can reproduce the complex behavior of the initial system. Because of the presence of a small perturbation parameter (i.e., the mass ratio between the BNES and the VdP oscillator), the APMD is governed by two different time scales. More precisely, it appears as a (3,1)-fast-slow system whose motion is constituted in a succession of slow and fast epochs. Founding a (3,1)-fast-slow APMD is interesting since that implies a more complex dynamics than in the case of a classic NES whose APMD is only (2,1)-fast-slow. A fast-slow analysis is finally conducted within the framework of the geometric singular perturbation theory. From the computation of the so-called critical manifold and the analytical expressions of the APMD fixed points, a global stability analysis is performed. This enables us to interpret a certain number of regimes observed on numerical simulations of the initial system.

**Keywords:** Passive vibration control, Bistable nonlinear energy sink, Self-sustained oscillations, Fast-slow analysis, Critical manifold

## C.1 INTRODUCTION

The well-known nonlinear energy sinks (NESs) are nonlinear vibration absorbers used for passive mitigation of unwanted oscillations caused by either free, external, parametric or self-excitations of a mechanical or acoustical system referred to as the primary structure. In most cases, an NES is a light mass (compared to the total mass of the primary structure) nonlinear attachment consisting of an essentially nonlinear spring (most of the time purely cubic) and a viscous linear damper. The dependence between the vibratory amplitude and the oscillating frequency of the NES (because of its strongly nonlinear nature) makes it able to resonate at any frequency. After tuning to the primary structure and absorbing and dissipating its energy, the NES can detune to prevent the energy from returning. This irreversible transfer of vibrational energy from the primary structure to the NES is called targeted energy transfer (TET) or energy pumping. In the seminal papers by Gendelman, Vakakis et al. [76, 203] the TET phenomenon is explained by the interaction between two nonlinear normal modes of vibration of the system producing a 1:1 resonance capture. Reviews of these concepts can be found in [206].

A few years later, the bistable NES (BNES), consisting of a small mass connected to a linear primary oscillator by a spring with both negative linear and nonlinear (cubic) stiffnesses, has been introduced by Gendelman and Lamarque [79]. In this seminar paper, various dynamical regimes were highlighted and conditions for efficient TET from the primary linear oscillator to the BNES (even though the term was not yet used there) were formulated. In particular, possibility of transient chaotic response was investigated. Other studies appeared much later [134, 171]. The negative linear and nonlinear stiffness components can be realized for example
through the geometric nonlinearity of the transverse displacement of pre-compressed linear springs [180] or
considering the transverse vibration of a thin viscously damped fixed-fixed beam with a small mass fixed at
its centre when buckled under axial constraint [103, 136, 170]. Most studies about BNESs consider either
mitigation of free vibrations or vibrations caused by external harmonic forcing. The first papers [79, 134, 171]
relate to mitigation of free vibrations. Romeo et al. [170] show that BNES can produce TET even at low
energies, in contrast to the classic cubic NESs. From an analytical method, Habib and Romeo [95] propose a
procedure to tune and optimize a BNES to mitigate free vibration of a two-degree-of-freedom (2-dof) primary
structure. Dekemele et al. [51] present an analytical study of the tuning and the performance of a BNES
again in the context of mitigation of free vibrations. Mattei et al. [136] and Iurasov and Mattei [103] study
both numerically and experimentally the mitigation of forced vibrations by a beam-based BNES. An adapted
semi-analytical method to predict the energy pumping time of a BNES under harmonic forcing is proposed
by Wu et al. [217]. The vibration reduction performance of a BNES and an NES is compared, by Wand et
al. [211] and the results show that the BNES performs over the full frequency band.

Using a cubic NES to mitigate limit cycle oscillations (LCOs) resulting from dynamic instabilities has been extensively studied in the literature. Mitigation of LCOs of the Van der Pol oscillator has been studied numerically by Lee et al. in their seminal paper [122] and then analytically by Gendelman and Bar [75]. A number of works focus on mitigation of flutter instabilities in aircraft wings. This problem has been first studied both numerically and experimentally, again by Lee. et al. [121, 123, 124]. The theoretical prediction of the observed operating regimes has been performed by means of multiple time scales approaches [78, 131] and improved more recently using the center manifold reduction technique [a2]. The possibility of mitigating self-sustained oscillations of a linear friction system having two unstable modes has been studied by Bergeot et al. [a9] by means of a sophisticated multiple time scales analysis allowing the understanding of the phenomena underlying the appearance of the many possible regimes of the system. Effect of stochastic forcing on the dynamic behavior of a self-sustained oscillator coupled to a nonlinear energy sink has been investigated in [a3]. On the contrary, the use of a BNES to mitigate self-sustained vibrations is poorly addressed in the literature. To the best of the authors' knowledge, the only paper dealing with this subject is by Franzini et al. [70]. It presents numerical and analytical studies for a BNES-based device used to simultaneously alleviate the unwanted vibrations caused by a galloping phenomenon and perform electro-mechanical energy conversion by the piezoelectric effect.

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. It is customary to study the dynamic behavior of the system in the neighborhood of this 1 :1
resonance capture [76, 203]. In general, the system is simplified by averaging it over a natural period of the
primary structure using the so-called Complexification-Averaging Method (CAM) [135, 206] or the Multiple
Scale/Harmonic Balance Method (MSHBM) [130, 228]. The resulting system, called here amplitude-phase
modulation dynamics<fnref n="1" /> (APMD), is then analyzed by means of singular perturbation techniques (by multiple
scales methods [144] or geometric singular perturbation theory [107]). The first key point of these analytical
processes is the partitioning of the APMD into two time scales (one slow and one fast <fnref n="2" />). In this representation
the time evolution of the APMD is thus described as a succession of slow and fast epochs which are ana-
lyzed independently. The second key point is the definition of the so-called critical manifold<fnref n="3" /> whose system
trajectories converge during slow epochs.

Some works mentioned above on the analytical study of energy pumping by a BNES follows this multiple time scales based methodology (or fast-slow analysis). However, in these studies, the CAM (or MSHBM) is used in the same way as if we were considering a classical NES, i.e., assuming that the BNES motion has the form of a zero-mean oscillating regime. This forgets one essential feature a bistable system motion, it can have a nonzero-mean. The resulting analysis can therefore only partially explain the dynamic behavior of the system, in particular for chaotic and intra-well motions of the BNES.

In the present paper, the MSHBM is used and modified (compared to its usual use) to take into account a
nonzero-mean motion. This allows us to obtain an APMD and a critical manifold that can explain numerous
regimes observed in numerical simulations of a Van der Pol oscillator (used as an archetypal self-sustained
oscillator) coupled to a BNES. Understanding the mechanisms underlying the responses observed in numerical
simulations of this simple model could enable us in the future to interpret the results of numerical simulations

<footnote n="1">1. Most of the time authors name it the *slow flow dynamics*. It is named here in the manner of Luongo and Zulli [130, 228] to avoid vocabulary ambiguities.</footnote>

<footnote n="2">2. Sometimes the terms “slow” and “super-slow” [59, 75] are used to leave the term “fast” for the oscillations on which the averaging was carried out and therefore to be coherent with the 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.</footnote>

<footnote n="3">3. The critical manifold is also called *slow invariant manifold* in the literature.</footnote>
<figure><img src="image_83.png" /><figcaption>FIGURE C.1. A Van der Pol oscillator coupled to a BNES.</figcaption></figure>

of more complex models and/or experimental data.

The paper is organized as follows. In Section C.2, the equations of motion of the full-order system under study - i.e. a Van der Pol oscillator coupled to a BNES - are obtained. Then, the dynamic behavior of a Van der Pol oscillator coupled to a BNES is investigated by means of numerical simulations. Using the MSHBM, the equations of the APMD are derived in Section D.3.1. The fast-slow analysis of the APMD is detailed in Section C.4. That allows us to perform a asymptotic global stability analysis in Section C.5, which enables us to understand the observed behavior using numerical simulations. Finally, concluding remarks and perspectives are formulated in Section C.6.

## C.2 EQUATIONS AND DYNAMIC BEHAVIOR OF THE MODEL

### C.2.1 Equations of motion of a Van der Pol oscillator coupled to a BNES

In this section we derive the equations of motion of the full-order system under study, i.e. a Van der Pol (VdP) oscillator - used as an archetypal self-sustained oscillator - coupled to an ungrounded bistable nonlinear energy sink (BNES). The BNES consists in a mass coupled to a linear damper and to a spring with a negative stiffness and a cubic stiffness. A schematic representation of this system is shown in Fig. E.1. The equations of motion of this system, with respect to the physical time $\hat{t}$, are as follows

$$ m_1 \frac{d^2 \hat{x}}{dt^2} - c_1 \frac{d\hat{x}}{dt} + c_1^{\text{NL}} \frac{d\hat{x}}{dt} \hat{x}^2 + k_1 \hat{x} + c_2 \left( \frac{d\hat{x}}{dt} - \frac{d\hat{y}}{dt} \right) - k_2 (\hat{x} - \hat{y}) + k_2^{\text{NL}} (\hat{x} - \hat{y})^3 = 0 \quad (C.1a) $$

$$ m_2 \frac{d^2 \hat{y}}{dt^2} + c_2 \left( \frac{d\hat{y}}{dt} - \frac{d\hat{x}}{dt} \right) - k_2(\hat{y} - \hat{x}) + k_2^{\text{NL}}(\hat{y} - \hat{x})^3 = 0 \quad (C.1b) $$

where $m_1$ and $m_2$ are the masses of the primary VdP oscillator and of the BNES, respectively. The parameter $k_1$ is the linear stiffness of the VdP oscillator and $-c_1$ (with $c_1 > 0$) and $c_1^{\text{NL}}$ characterize its linear negative and nonlinear dampings, respectively. The BNES is characterized by its linear damping coefficient $c_2$, negative linear stiffness $-k_2$ (with $k_2 > 0$) and nonlinear stiffness $k_2^{\text{NL}}$.

Equation (G.8) is rescaled leading to the following dimensionless system of differential equations

$$ \ddot{x} - \epsilon\sigma\dot{x} + \epsilon\lambda\dot{x}x^2 + x + \epsilon\mu(\dot{x} - \dot{y}) - \epsilon\beta(x - y) + \epsilon\alpha(x - y)^3 = 0 \quad (C.2a) $$

$$ \epsilon\ddot{y} + \epsilon\mu(\dot{y} - \dot{x}) - \epsilon\beta(y - x) + \epsilon\alpha(y - x)^3 = 0 \quad (C.2b) $$

where $x = \hat{x}/L$, $y = \hat{y}/L$ (with $L$ a characteristic length of the motion), $\omega_1 = \sqrt{k_1/m_1}$, $t = \omega_1\hat{t}$, $\{\} = d\{\}/dt$, $\epsilon = m_2/m_1$ is the mass ratio between the NES and the VdP oscillator, $\sigma = c_1/(m_2\omega_1)$, $\lambda = c_1^{\text{NL}}L^2/(m_2\omega_1)$, $\mu = c_2/(m_2\omega_1)$, $\beta = k_2/(m_2\omega_1^2)$ and $\alpha = k_2^{\text{NL}}L^2/(m_2\omega_1^2)$.

The static solutions of Eq. (E.3) are : the trivial solution

$$ p_0^{\text{st}} = (x_0^{\text{st}}, y_0^{\text{st}}) = (0, 0) \quad (C.3) $$
which is unstable and the two following non trivial solutions

$$p_1^{\text{st}} = (x_1^{\text{st}}, y_1^{\text{st}}) = \left(0, \sqrt{\frac{\beta}{\alpha}}\right) \quad \text{and} \quad p_2^{\text{st}} = (x_2^{\text{st}}, y_2^{\text{st}}) = \left(0, -\sqrt{\frac{\beta}{\alpha}}\right). \qquad (C.4)$$

which are stable if there are no self-sustained oscillations. For a classic cubic NES, i.e., $\beta = 0$, the only equilibrium position is $p_0^{\text{st}}$ which is stable (again if there are no self-sustained oscillations).

Finally, for comparison purposes, we give here an approximated solution of the VdP oscillator uncoupled from the NES whose equation of motion is $\ddot{x} - \epsilon\sigma\dot{x} + \epsilon\lambda\dot{x}x^2 + x = 0$. This approximated solution is obtained using the Krylov-Bogoliubov method of averaging (see e.g. [143]) as

$$x(t) = \frac{2\sqrt{\sigma}e^{\frac{\sigma t\epsilon}{2}}}{\sqrt{\frac{4\sigma}{r_0^2} + \lambda(e^{\sigma t\epsilon} - 1)}} \cos(t + \varphi_0) \qquad (C.5)$$

where $r(t)$ is the amplitude of $x(t)$, $r_0 = r(t = 0)$ and $\varphi_0$ is the phase at the origin of time. Details on obtaining Eq. (C.5) are given in Appendix C.7. For $\lambda > 0$, the amplitude $r(t)$ in Eq. (C.65) tends to zero if $\sigma < 0$ and to $2\sqrt{\frac{\sigma}{\lambda}}$ if $\sigma > 0$. In other words, for the original VdP oscillator (C.58), if $\sigma < 0$ the trivial solution $(x, \dot{x}) = (0, 0)$ is stable and becomes unstable when $\sigma > 0$. The periodic solution

$$(x, \dot{x}) = \left( 2\sqrt{\frac{\sigma}{\lambda}} \cos(t + \varphi_0), -2\sqrt{\frac{\sigma}{\lambda}} \sin(t + \varphi_0) \right) \qquad (C.6)$$

exists and it is stable for $\sigma > 0$.

### C.2.2 Response regimes of a Van der Pol oscillator coupled to a BNES

In this section, the dynamic behavior of a Van der Pol oscillator coupled to a BNES is presented by means of direct numerical integration of Eq. (E.3).

We first recall briefly the response regimes observed when an NES (i.e., with no linear stiffness) is used instead of a BNES. As it is widely discussed in the literature (see e.g. [75, 78]), in this case, four scenarios are observed. In previous works by the authors [a7, a8] these responses are sorted into two categories. In the first category, the NES acts, resulting in three responses called *harmless situations*. These three harmless regimes are: (i) the *complete suppression* in which the trivial equilibrium is stabilized because of the additional linear part of the NES including mass and damping; (ii) the *mitigation through periodic responses (PRs)* with an amplitude smaller than the periodic response undergone by the VdP oscillator alone and (iii) the *mitigation through a Strongly Modulated Response (SMR)* corresponding to a quasi-periodic (amplitude and phase modulated) regime again with a maximum amplitude smaller than the periodic response undergone by the VdP oscillator alone. In the second category, the NES does not act, resulting in a *no mitigation response* called *harmful situation*. In this case a periodic regime with an amplitude close to what is observed when the VdP oscillator is not coupled to the NES. The specific value of the bifurcation parameter $\sigma$ at which the system switches from a harmless regime to a harmful regime is called the *mitigation limit* of the NES, denoted as $\sigma_{ml}$. In the case of an NES the mitigation limit is defined considering a set of initial conditions (for the slow flow) as a small perturbation of the trivial equilibrium solution, i.e., $(x_0^{\text{st}}, y_0^{\text{st}}) = (0, 0)$.

The specific features of a BNES compared with an NES are presented first in Fig. C.2 in which, as a function of $\sigma$, the steady-state maximum amplitude of the VdP oscillator alone is compared to the steady-state maximum amplitude of the VdP oscillator coupled to an NES on the one hand and to a BNES on the other. The following set of parameters is used:

$$\text{VdP parameter: } \lambda = 0.5 \qquad (C.7a)$$

$$\text{BNES parameters: } \mu = 0.06, \quad \beta = 0.575, \quad \alpha = 2 \qquad (C.7b)$$

$$\text{Mass ratio: } \epsilon = 0.0025 \qquad (C.7c)$$

The steady-state maximum amplitude of the VdP oscillator alone is defined as follows

$$A_x^{\text{wo}} = \frac{\max[(x(0.9T), x(T))] - \min[(x(0.9T), x(T))]}{2} \qquad (C.8)$$

where $x(t)$ is the time series of the variable $x$ computed by means of the numerical integration of Eq. (C.58) from $t = 0$ to $t = T$ and using the parameters (C.7a). In a similar way, the steady-state maximum amplitude
<figure><img src="image_84.png" /><figcaption>FIGURE C.2. Steady-state maximum amplitude of the VdP oscillator alone $A_x^{wo}$ (green) compared to the steady-state maximum amplitudes $A_x^{wNES}$ and $A_x^{wBNES}$ of the VdP oscillator coupled to an NES (red) and to a BNES (blue), respectively, as functions of the bifurcation parameter $\sigma$. The theoretical amplitude of the VdP alone (see Eq. (C.6)) is also displayed (black dashed). Parameters: (C.7a) are used to compute $A_x^{wo}$ from Eq. (C.58); (C.7a) with $\mu = 0.4, \beta = 0$ are used to compute $A_x^{wNES}$ from Eq. (E.3) and (C.7) are used to compute $A_x^{wBNES}$ also from Eq. (E.3). The value of $\epsilon$ is given by Eq. (C.7c).</figcaption></figure>

of the VdP oscillator coupled to an NES is

$$A_x^{wNES} = \frac{\max[(x(0.9T), x(T))] - \min[(x(0.9T), x(T))]}{2} \qquad (C.9)$$

where $x(t)$ is the time series of the variable $x$ computed by means of the numerical integration of Eq. (E.3) with $\mu = 0.4, \beta = 0$ and $\alpha = 2$ (case of an NES with good properties of mitigation) and using the VdP parameters (C.7b). Finally, the steady-state maximum amplitude of the VdP oscillator coupled to a BNES is

$$A_x^{wBNES} = \frac{\max[(x(0.9T), x(T))] - \min[(x(0.9T), x(T))]}{2} \qquad (C.10)$$

where $x(t)$ is the time series of the variable $u$ computed by means of the numerical integration of Eq. (E.3) using the parameters (C.7) (case of an NES with good properties of mitigation).

The theoretical amplitude of the VdP alone (see Eq. (C.6)) is also displayed in Fig. C.2. We can first see that the mitigation limit of the BNES ($\sigma_{ml} \approx 8.03$) is higher than that of the NES ($\sigma_{ml} \approx 2.03$). Secondly, in cases of attenuation (PRs or SMRs, the difference between these two regimes is not visible in the figure: from $\sigma \approx 0.39$ to $\sigma \approx 2.03$ for the NES and from $\sigma \approx 1.87$ to $\sigma \approx 8.03$ for the BNES) the maximum amplitude $A_x^{wBNES}$ is almost 10 times smaller than $A_x^{wNES}$. This last observation is of course interesting from the point of view of vibration attenuation, but we shall see that it can have a negative impact on the robustness of the BNES.

How large is the perturbation does not really matter in the case a classic purely cubic NES. This is not the case for a BNES. To show this, Fig. C.3 plots the mitigation limit $\sigma_{ml}$ as a function of $\beta$ and $\mu$ in the case of an NES and a BNES for $\lambda$ given by (C.7a), $\alpha = 2$ and $\epsilon$ given by Eq. (C.7c). The initial conditions are $(x_0, y_0) = p_0^\text{st} + (\text{per}, \text{per})$ for the NES and $(x_0, y_0) = p_1^\text{st} + (\text{per}, \text{per})$ for the BNES, respectively. The perturbation value is $\text{per} = 0.001, 0.01$ and $0.1$ in Figs. C.3(a) to C.3(c), respectively. The mitigation limit $\sigma_{ml}$ is computed by locating the jumps on amplitude profiles such as those shown in Fig. C.2. In the latter, results have been obtained in the same way with $\text{per} = 0.001$.

As mentioned above, the mitigation limit of the NES (in blue color in Fig. C.3) is not affected by the change in initial conditions and of course it does not depend on $\beta$. On the contrary, the mitigation limit of the BNES (in orange color in Fig. C.3) it is greatly affected by the change in initial conditions. Indeed, a very high maximum is present when $\text{per}$ is small. This maximum decreases until it disappears when $\text{per}$ increases. That shows that the BNES can be very efficient but this is not robust since a little perturbation can switch the system from harmless to harmful situations. However, even in the case of $\text{per} = 0.1$ (see Fig. C.3(c)) it is possible to find a set of parameters for which the BNES performs better than the NES.

The overall results of Fig C.3, in terms of mitigation limit, do not highlight the complex dynamic behavior of the BNES. As in the case of a classic NES, when $\sigma > \sigma_{ml}$ the system (E.3) undergoes a periodic regime
<figure><img src="image_85.png" /><figcaption>**FIGURE C.3.** Mitigation limit $\sigma_{ml}$ as a function of $\beta$ and $\mu$ in the case of an NES (blue surface) and a BNES (orange surface) for $\lambda$ given by (C.7a), $\alpha = 2$ and $\epsilon$ given by Eq. (C.7c). The initial conditions are $(x_0, y_0) = p_0^{\text{st}} + (\text{per}, \text{per})$ for the NES and $(x_0, y_0) = p_1^{\text{st}} + (\text{per}, \text{per})$ for the BNES, respectively. The perturbation parameter is 0.001, 0.01 and 0.1 in (a) to (c), respectively.</figcaption></figure>

with an amplitude for the variable $x$ close to that of the VdP oscillator alone. However, when $\sigma < \sigma_{ml}$ a multitude of harmless regimes can be observed, whereas there are only three responses in the case of an NES. In Fig. C.4, seven of these regimes are presented. The figure shows the time series $y(t)$ obtained from the numerical simulation of Eq. (E.3) and only steady states are shown. The equilibrium positions $y_1^{\text{st}}$ and $y_2^{\text{st}}$ (see Eq. (C.4)) are also depicted. The parameters used are $\lambda = 0.5$, $\alpha = 2$, $\epsilon = 0.0025$ (as in Fig. C.3) and those given in Tab. C.1. The initial conditions are the same as in Fig. C.3(a), i.e., with $\text{per} = 0.001$. From Fig. C.4(a) to Fig. C.4(g) one has:

* Fig. C.4(a). The stabilization of $p_2^{\text{st}}$.

* Fig. C.4(b). A PR centered on $y = 0$ (PR1).

* Fig. C.4(c). An SMR centered on $y = 0$ (SMR1).

* Fig. C.4(d). An intra-well (i.e., centered on $y = y_1^{\text{st}}$) PR (PR2).

* Fig. C.4(e). An intra-well SMR (SMR2).

* Fig. C.4(f). An intermittent chaotic regime of type 1 (CR1), i.e., with a succession of non-chaotic oscillations centered on $y = 0$, chaotic inter-well oscillations and non-chaotic oscillations centered on $y = y_1^{\text{st}}$ or $y = y_2^{\text{st}}$. A zoom of Fig. C.4(f) is shown in Fig. C.5(a). On it, we can see first non-chaotic oscillations around $y = y_2^{\text{st}}$, then a chaotic motion, then non-chaotic oscillations around zero and again non-chaotic oscillations around $y = y_2^{\text{st}}$. Concerning these oscillations around $y = y_2^{\text{st}}$, we first see amplitude-modulated oscillations that abruptly diminish to almost zero; after that we see oscillations with exponential growth before returning to a chaotic movement. The phase in which the oscillations around $y = y_2^{\text{st}}$ (or also $y = y_1^{\text{st}}$) are amplitude-modulated only rarely appears.

* Fig. C.4(g). An intermittent chaotic regime of type 2 (CR2), i.e., with a succession of non-chaotic oscillations centered on $y = 0$ and chaotic inter-well oscillations. This appears more clearly in Fig. C.5(b) showing a zoom of Fig. C.4(g).

For regimes of Fig. C.4(d) and (e), the motions observed are centered on $y_1^{\text{st}}$ but similar regimes centered on $y_2^{\text{st}}$ can also be observed depending on the chosen initial conditions. Regimes (a) to (c) are also observed with a classic cubic NES (of course in this case the stabilization concerns the trivial solution $y = y_0^{\text{st}}$) while the regimes (d) to (g) are specific to BNES.

The maximum of $\sigma_{ml}$ observed in Fig. C.3(a) corresponds to a stabilization of the equilibrium $y_1^{\text{st}}$ or $y_2^{\text{st}}$ and for the other high efficient situations, the regimes observed just before $\sigma_{ml}$ are the SMR2s or less often PR2s (for even smaller values of the perturbation parameter $\text{per}$, these two regimes with very small amplitudes can correspond to the highest values of $\sigma_{ml}$). Responses of type SMR2s are observed in the most efficient situations in the cases of Figs. C.3(b) and C.3(c). These observations will be interpreted after the fast-slow analysis carried out in Section C.4. Before that, the equation of the so-called *amplitude-phase modulation* dynamics are derived in the next section.
<figure><img src="image_86.png" /><figcaption>**FIGURE C.4.** Time series $y(t)$ (solid blue line) obtained from the numerical simulation of Eq. (E.3), only steady states are shown. Equilibrium positions $y_1^{st}$ (dashed green line) and $y_2^{st}$ (dashed red line). From (a) to (g) one has: (a) the stabilization of $p_2^{st}$; (b) a PR centered on $y=0$ (PR1); (c) an SMR centered on $y=0$ (SMR1); (d) a PR centered on $y=y_2^{st}$ (PR2); (e) an SMR centered on $y=y_2^{st}$ (SMR2); (f) an intermittent chaotic regime of type 1 (CR1), i.e., with a succession of oscillations centered on $y=0$ and chaotic motions; and (g) an intermittent chaotic regime of type 2 (CR2), i.e., with succession of oscillations centered on $y=0$, chaotic motions and oscillations centered on $y=y_1^{st}$ or $y=y_2^{st}$. The parameters used are $\lambda=0.5$, $\alpha=2$, $\epsilon=0.0025$ and those given in Tab. C.1. The initial conditions are the same as in Fig. C.3(a), i.e., with $per=0.001$.</figcaption></figure>

<table><caption>**TABLEAU C.1.** Parameters used for numerical simulations shown in Fig. C.4. The corresponding values of $\sigma_{ml}$ are also indicated.</caption><thead><tr><th></th><th>β</th><th>μ</th><th>σ</th><th>σ<sub>ml</sub></th></tr></thead><tbody><tr><td>(a)</td><td>0.52</td><td>0.19</td><td>3.9</td><td>5.2</td></tr><tr><td>(b)</td><td>0.12</td><td>0.55</td><td>0.78</td><td>1.73</td></tr><tr><td>(c)</td><td>0.12</td><td>0.31</td><td>1.22</td><td>1.63</td></tr><tr><td>(d)</td><td>0.62</td><td>0.19</td><td>3.75</td><td>5</td></tr><tr><td>(e)</td><td>0.62</td><td>0.07</td><td>3.96</td><td>5.28</td></tr><tr><td>(f)</td><td>0.92</td><td>0.07</td><td>1.43</td><td>1.9</td></tr><tr><td>(g)</td><td>1</td><td>0.19</td><td>1.74</td><td>2.33</td></tr></tbody></table>
<figure><img src="image_87.png" /><figcaption>FIGURE C.5. Zooms of (a) Figs. C.4(f) and (b) C.4(g).</figcaption></figure>

C.3 EQUATIONS OF THE AMPLITUDE-PHASE MODULATION DYNAMICS

TET is due to the interaction between two nonlinear modes of the coupled structure [76, 203]. This phenomenon, called a 1 :1 resonance capture, occurs at a frequency close to the natural frequency of the primary structure, here the VdP oscillator. It is customary to study the dynamic behavior of the system in the neighborhood of this 1 :1 resonance capture. In general, the system is simplified by averaging it over a natural period of the primary structure using the so-called Complexification-Averaging Method (CAM) [135, 206]. The resulting averaged dynamics is called amplitude-phase modulation dynamics (APMD). The CAM is not consistent here due the bistable nature of the BNES, i.e., the fact that its motion can have a nonzero-mean. This is the reason why the Multiple Scale/Harmonic Balance Method (MSHBM) [130, 228] is preferred because it can easily and rigorously be modified (compared to its usual use) to take into account a nonzero-mean motion.

First, using the following change of variable

$$
u = x + \epsilon \quad \text{and} \quad v = x - y \tag{C.11}
$$

the equations of motion (E.3) become

$$
\ddot{u} + u + \epsilon [-u - \sigma \dot{u} + \lambda \dot{u} u^2 + v] = 0 \quad (C.12a)
$$

$$
\ddot{v} + \mu \dot{v} - \beta v + \alpha v^3 + u = 0 \quad (C.12b)
$$

Assuming a small mass ratio between the NES and the Van der Pol oscillator (i.e., 0 < ϵ ≪ 1), only terms of order equal to 1 or less in ϵ have been kept in Eq. (C.12a). In Eq. (C.12b) keeping only the terms of order 0 is sufficient to correctly reproduce the essential characteristics of the dynamic behavior of the system.

Following the Multiple Scale Method (MSM) [143] the independent time variables $T_0 = t$, $T_1 = \epsilon t$, $T_2 = \epsilon^2 t$, ... are introduced. The time derivative rules are therefore $\frac{d}{dt} = D_0 + \epsilon D_1 + \epsilon^2 D_2 + ...$ and $\frac{d^2}{dt^2} = D_0^2 + 2\epsilon D_0 D_1 + \epsilon^2 (D_0 D_2 + D_1^2) + ...$ where $D_i = \frac{\partial}{\partial T_i}$ ($i = 0, 1, 2, ...$). The variables $u$ and $v$ are expanded in power series of $\epsilon$ as

$$
u = u_0 + \epsilon u_1 + \epsilon^2 u_2 + \dots \tag{C.13a}
$$

$$
v = v_0 + \epsilon v_1 + \epsilon^2 v_2 + \dots \tag{C.13b}
$$

In the sequel, only the first two time scales $T_0$ and $T_1$ are retained. Substituting Eq. (C.13) and the time derivative rules into Eq. (E.5) and equating coefficients of same powers of $\epsilon$, we obtain

$$
\text{Order } \epsilon^0: \quad D_0^2 u_0 + u_0 = 0 \tag{C.14a}
$$

$$
D_0^2 v_0 + \mu D_0 v_0 - \beta v_0 + \alpha v_0^3 = 0. \quad (C.14b)
$$

$$
\text{Order } \epsilon^1: \quad D_0^2 u_1 + u_1 = -\lambda u_0^2 D_0 u_0 + \sigma D_0 u_0 - 2 D_0 D_1 u_0 + u_0 - v_0 \quad (C.15a)
$$
$$D_0^2 v_1 + \mu D_0 v_1 - \beta v_1 + \alpha v_1^3 = - \mu D_1 v_0 - 2 D_0 D_1 v_0 + u_0 - u_1 + \alpha v_1^3 - 3 \alpha v_0^2 v_1. \quad (C.15b)$$

The harmonic solution of Eq. (C.14a) is written on the following complex form

$$u_0(T_0, T_1) = \frac{A(T_1)e^{iT_0} + A^*(T_1)e^{-iT_0}}{2} \quad (C.16)$$

where $A^*$ is the complex conjugate of the complex amplitude $A$.

Equation (C.14b) is the equation of motion of the BNES alone and has no analytical solution. Then the Harmonic Balance Method (HBM) is used to find an approximated solution as the following two-terms Fourier series

$$v_0(T_0, T_1) = b_0(T_1) + \frac{C_0(T_1)e^{iT_0} + C_0^*(T_1)e^{-iT_0}}{2} \quad (C.17)$$

where the first term in the right-hand side is real and used to consider the bistable nature of the BNES, i.e., the fact that $v_0$ can have nonzero-mean. The second term indicates the simplification - due to the assumption of a 1:1 resonance capture - of the oscillating motion of $v_0$ as a simple harmonic motion.

Substituting Eq. (C.17) into Eq. (C.14b) and equalizing the coefficients of $e^0$ (mean terms) and $e^{iT_0}$ (first harmonic terms) the following complex algebraic equations are obtained,

$$-\beta b_0 + \frac{3}{2} \alpha C_0 |b_0|^2 + \alpha b_0^3 = 0 \quad (C.18a)$$

$$4 (A + C_0 (-\beta + 3\alpha b_0^2 + i\mu - 1)) + 3\alpha C_0 |C_0|^2 = 0, \quad (C.18b)$$

respectively.

Let us now deal with Eq. (C.15). The first step is to substitute Eqs. (C.16) and (C.17) into Eq. (C.15a). Secular terms, i.e., terms in factor of $e^{iT_0}$ and $e^{-iT_0}$ in the right-hand side of the resulting equation, create diverging solutions and must be eliminated. These terms in factor of $e^{iT_0}$ and $e^{-iT_0}$ correspond to two equivalent (i.e., complex conjugate) differential equations with respect $D_1 A$ and $D_1 A^*$, respectively, and called solvability conditions. By arbitrarily selecting the first condition, one obtains the following equation

$$D_1 A = \frac{1}{8} (4A(\sigma - i) + 4iC_0 - \lambda A |A|^2). \quad (C.19)$$

The particular non-diverging solution of Eq. (C.15a) is then obtained, it is

$$u_1(T_0, T_1) = -b_0(T_1) - \frac{1}{64} i \lambda (A(T_1)^3 e^{3iT_0} - A^*(T_1)^3 e^{-3iT_0}). \quad (C.20)$$

The second step is to find $v_1(T_1, T_2)$ and again the HBM is used assuming the following form

$$v_1(T_0, T_1) = b_1(T_1) + \frac{C_1(T_1)e^{iT_0} + C_1^*(T_1)e^{-iT_0}}{2}. \quad (C.21)$$

Substituting Eqs. (C.16), (C.17), (C.20) and (C.21) into Eq. (C.15b) and again equating the coefficients of $e^0$ and $e^{iT_0}$ the following differential equations are obtained

$$D_1 b_0 = \frac{1}{2\mu} \left( b_1 \left( 2\beta - 3\alpha C_0 \overline{C_0} \right) + b_0 \left( 2 - 3\alpha \left( C_0 \overline{C_1} + C_1 \overline{C_0} \right) \right) - 6\alpha b_1 b_0^2 \right), \quad (C.22)$$

and

$$D_1 C_0 = \frac{1}{4(\mu + 2i)} \left( 2C_1 \left( -3\alpha C_0 \overline{C_0} + 2\beta - 6\alpha b_0^2 - 2i\mu + 2 \right) - 3\alpha C_0 \left( C_0 \overline{C_1} + 8b_0 b_1 \right) \right), \quad (C.23)$$

respectively.

The MSM procedure is stopped at order $\epsilon^1$. The last step of the MSHBM is then to come back to the physical time $t$. To achieve that we use the rules $b = b_0 + \epsilon b_1$ and $C = C_0 + \epsilon C_1$ and the fact that the variables $A$, $B$ and $C$ do not depend on the fast time $T_0$, therefore for them we have $\frac{d}{dt} = \epsilon D_1$. First, from Eq. (C.19) we obtain

$$\dot{A} = \frac{\epsilon}{8} (4A(\sigma - i) + 4iC - \lambda A |A|^2). \quad (C.24)$$

Then combining Eq. (C.18a) with Eq. (C.22) and Eq. (C.18b) with Eq. (C.23) and keeping only terms of order equal to 1 or less in $\epsilon$ we have

$$\dot{b} = \frac{b(2\beta - 2\alpha b^2 - 3\alpha |C|^2)}{2\mu} + \epsilon \frac{b}{\mu}, \quad (C.25)$$
and

$$ \dot{C} = \frac{C(-3\alpha(|C|^2 + 4b^2) + 4(\beta - i\mu + 1)) - 4b}{4(\mu + 2i)}, \qquad (C.26) $$

respectively.

Equations (C.24) to (C.26) constitute the complex form of the APMD. Finally, substituting A and C in the complex APMD by their polar coordinates defined as $A = ae^{j\theta}$ and $C = ce^{j\phi}$, new equations of motion for the real amplitudes a, b and c and the phase difference $\delta = \theta - \phi$ is obtained as

$$ \dot{a} = \epsilon f(a, c, \delta) \qquad (C.27a) $$

$$ \dot{b} = g_1(b, c, \epsilon) \qquad (C.27b) $$

$$ \dot{c} = g_2(a, b, c, \delta) \qquad (C.27c) $$

$$ \dot{\delta} = g_3(a, b, c, \delta, \epsilon) \qquad (C.27d) $$

with

$$ f(a, c, \delta) = \frac{1}{8}(4a\sigma - a^3\lambda + 4c\sin\delta) \qquad (C.28a) $$

$$ g_1(b, c, \epsilon) = \frac{b(2\beta - 2\alpha b^2 - 3\alpha c^2)}{2\mu} + \epsilon \frac{b}{\mu} \qquad (C.28b) $$

$$ g_2(a, b, c, \delta) = -\frac{4a\mu \cos\delta + 8a \sin\delta + c\mu(-4\beta + 3\alpha(4b^2 + c^2) + 4)}{4(\mu^2 + 4)} \qquad (C.28c) $$

$$ g_3(a, b, c, \delta, \epsilon) = -\frac{4a \cos \delta - 2a\mu \sin \delta + 3\alpha c (4b^2 + c^2) - 2c (2\beta + \mu^2 + 2)}{2c (\mu^2 + 4)} - \frac{\epsilon a - c \cos \delta}{2a}. \qquad (C.28d) $$

Equation (C.27) constitutes the real form of the APMD.

Due to the presence of the small parameter $\epsilon$, the APMD is governed by two different time scales : the fast time scale $t$ and the slow time scale $\tau = \epsilon t$ (note that $t$ and $\tau$ were denoted above $T_0$ and $T_1$, respectively, in the context of the MSM). More precisely, within the framework of the geometric singular perturbation theory [107], Eq. (C.27) appears as a (3, 1)-fast-slow system where $b, c$ and $\delta$ are the fast variables and $a$ the slow variable. Periodic steady-state regimes of the original mechanical system (E.3) correspond to non trivial fixed points of the APMD whereas periodic responses of the latter characterize amplitude and phase modulated responses of the original dynamics.

# C.4 FAST-SLOW ANALYSIS OF THE AMPLITUDE-PHASE MODULATION DYNAMICS

## C.4.1 The critical manifold

A key mathematical tool for the description of a fast-slow system such as Eq. (C.27) is its critical manifold. The latter is defined in this section. First, the APMD (C.27) is written with respect to the slow time $\tau$ as follows

$$ a' = f(a, c, \delta) \qquad (C.29a) $$

$$ \epsilon b' = g_1(b, c, \epsilon) \qquad (C.29b) $$

$$ \epsilon c' = g_2(a, b, c, \delta) \qquad (C.29c) $$

$$ \epsilon \delta' = g_3(a, b, c, \delta, \epsilon) \qquad (C.29d) $$

where $(\cdot)'$ denotes the derivative with respect to the slow time $\tau$.

Both in the forms of (C.27) and (C.29), the APMD is a (3, 1)-fast-slow system whose the time evolution is described by a succession of slow and fast epochs. Because $0 < \epsilon \ll 1$, the behavior of the APMD during each of these epochs can be approximated by that of the slow and fast subsystems defined as follows : considering $\epsilon = 0$ respectively in Eqs. (C.29) and (C.27) yields the slow subsystem

$$ a' = f(a, c, \delta) \qquad (C.30a) $$

$$ 0 = g_1(a, b, c, \delta, 0) \qquad (C.30b) $$

$$ 0 = g_2(a, b, c, \delta) \qquad (C.30c) $$
$$0 = g_3(a, b, c, \delta, 0) \qquad (C.30d)$$

which is a differential-algebraic equation, and the fast subsystem

$$\dot{a} = 0 \qquad (C.31a)$$

$$\dot{b} = g_1(a, b, c, \delta, 0) \qquad (C.31b)$$

$$\dot{c} = g_2(a, b, c, \delta) \qquad (C.31c)$$

$$\dot{\delta} = g_3(a, b, c, \delta, 0). \qquad (C.31d)$$

The critical manifold of the APMD (C.27) is the solution of the algebraic part of the slow sub system (C.30) and it is expressed as follows

$$\mathcal{M}_0 = \left\{ (a, b, c, \delta) \in \mathbb{R}^{+3} \times [-\pi, \pi] \mid g_1(b, c, 0) = 0, g_2(a, b, c, \delta) = 0 \text{ and } g_3(a, b, c, \delta, 0) = 0 \right\}. \quad (C.32)$$

The functions $g_1(b, c, 0) = 0$, $g_2(a, b, c, \delta) = 0$ and $g_3(a, b, c, \delta, 0) = 0$ are written as follows

$$g_1(a, b, c, \delta, 0) = b(b + H_{01}(c))(b - H_{01}(c)) \qquad (C.33a)$$

$$g_2(a, b, c, \delta) = H_{02}(b, c) + a(k_1 \cos \delta + k_2 \sin \delta) \qquad (C.33b)$$

$$g_3(a, b, c, \delta, 0) = H_{03}(b, c) + \frac{a}{c}(k_2 \cos \delta + k_1 \sin \delta). \qquad (C.33c)$$

where the function $H_{01}$ is given below in Eq. (C.38b) and the functions $H_{02}$ and $H_{03}$ and the parameters $k_1$ and $k_2$ can be easily deduced from Eqs. (C.28c) and (C.28d).

The critical manifold can be expressed as a one-dimensional manifold evolving in the four-dimensional phase space of the APMD. Indeed, solving first Eqs. (C.30c) and (C.30d) yield

$$\begin{aligned} \cos \delta &= -\frac{k_1 H_{02}(b, c) + k_2 c H_{03}(b, c)}{a(k_1^2 + k_2^2)} \\ &= -\frac{12\alpha b^2 c + 3\alpha c^3 - 4\beta c - 4c}{4a} \end{aligned} \qquad (C.34)$$

$$\begin{aligned} \sin \delta &= -\frac{k_2 H_{02}(b, c) - k_1 c H_{03}(b, c)}{a(k_1^2 + k_2^2)} \\ &= -\frac{c\mu}{a}. \end{aligned} \qquad (C.35)$$

Afterwards, combining Eqs. (C.34) and (C.35) with $a > 0$ leads to the following amplitude and phase equations

$$\begin{aligned} a &= \frac{\sqrt{c^2 H_{03}(b, c)^2 + H_{02}(b, c)^2}}{\sqrt{k_1^2 + k_2^2}} \\ &= \frac{1}{4} c \sqrt{(3\alpha(4b^2 + c^2) - 4(\beta + 1))^2 + 16\mu^2} \end{aligned} \qquad (C.36)$$

$$\begin{aligned} \tan \delta &= \frac{k_2 H_{02}(b, c) - k_1 c H_{03}(b, c)}{k_1 H_{02}(b, c) + k_2 c H_{03}(b, c)} \\ &= \frac{4\mu}{3\alpha(4b^2 + c^2) - 4(\beta + 1)}. \end{aligned} \qquad (C.37)$$

Then solving Eq. (C.30b) we obtain the three possible expressions of $b$ (which reflect the bistability of the BNES) as functions of $c$

$$b = 0 \qquad (C.38a)$$

$$b = \pm H_{01}(c) = \pm \sqrt{\frac{2\beta - 3\alpha c^2}{2\alpha}} \qquad (C.38b)$$

Equation (C.38b) shows that we must have $c < c^1$, with

$$c^1 = \sqrt{\frac{2\beta}{3\alpha}}, \qquad (C.39)$$
for b to be real. Note that by definition, substituting of Eqs. (C.34) and (C.35) into Eqs. (C.33b) and (C.33c), reduces the fast subsystem (C.31) to the following form

$$ \dot{a} = 0 \qquad (C.40a) $$

$$ \dot{b} = \frac{b(2\beta - 2\alpha b^2 - 3\alpha c^2)}{2\mu} = G(b, c) \qquad (C.40b) $$

$$ \dot{c} = 0 \qquad (C.40c) $$

$$ \dot{\delta} = 0. \qquad (C.40d) $$

It is then easy to show that (C.40b) undergoes a pitchfork bifurcation at ($b=0, c=c^I$), i.e., $G(0, c^I) = 0$, $\frac{\partial}{\partial b}G(0, c^I) = 0$, $\frac{\partial}{\partial c}G(0, c^I) = 0$, $\frac{\partial^2}{\partial b^2}G(0, c^I) = 0$, $\frac{\partial^2}{\partial b \partial c}G(0, c^I) \neq 0$ and $\frac{\partial^3}{\partial b^3}G(0, c^I) \neq 0$ (see e.g. [214] Chap. 20, for the definition of a pitchfork bifurcation).

Finally, combining Eq. (C.36) with each of the equations of (C.38) we obtain the two amplitudes branches of $M_0$ (only two because b appears squared in Eq. (C.36))

$$ a = H_1(c) = \frac{1}{4} c \sqrt{(3\alpha c^2 - 4(\beta + 1))^2 + 16\mu^2} \qquad (C.41a) $$

$$ a = H_2(c) = \frac{1}{4} c \sqrt{(-8\beta + 15\alpha c^2 + 4)^2 + 16\mu^2}. \qquad (C.41b) $$

In the same way, combining Eq. (C.37) with each of the equations of (C.38) we obtain

$$ \sin \delta = G_1^s(c) = -\mu_a^c, \qquad \cos \delta = G_1^c(c) = \frac{c(4\beta - 3\alpha c^2 + 4)}{4a}, \qquad (C.42a) $$

$$ \sin \delta = G_2^s(c) = -\mu_a^c, \qquad \cos \delta = G_2^c(c) = \frac{c(-8\beta + 15\alpha c^2 + 4)}{4a}. \qquad (C.42b) $$

The critical manifold can be therefore expressed as $M_0 = M_{01} \cup M_{02}$ with

$$ M_{01} = \left\{ (a, b, c, \delta) \in \mathbb{R}^{+3} \times [-\pi, \pi] \mid a = H_1(c), b = 0, \sin \delta = G_1^s(c), \cos \delta = G_1^c(c) \right\}, \qquad (C.43a) $$

$$ M_{02} = \left\{ (a, b, c, \delta) \in \mathbb{R}^{+2} \times [0, \sqrt{\frac{2\beta}{3\alpha}}] \times [-\pi, \pi] \mid a = H_2(c), b = \pm H_{01}(c), \sin \delta = G_2^s(c), \cos \delta = G_2^c(c) \right\}. \qquad (C.43b) $$

Note that the two branches intersect when $H_{01}(c) = 0$, i.e., for $c = c^I$ (see Eq. (C.39)).

Exploiting the polynomial properties of $H_1(c)$ and $H_2(c)$, the local extrema of these functions, i.e, the solutions of $\frac{dH_1}{dc} = 0$ and $\frac{dH_2}{dc} = 0$, respectively, can be computed. The local extrema of $H_1(c)$, denoted as $c_1^{LF}$ and $c_1^{RF}$, are determined (the explicit expressions are given by Eq. (C.66) in Appendix C.8) and exist if the following relation holds

$$ \mu < \frac{\beta + 1}{\sqrt{3}}. \qquad (C.44) $$

If the condition (C.44) is not satisfied, the function $H_1(c)$ no longer has local extrema. In the rest of the paper, one considers that the condition (C.44) always holds. The superscripts ()$^{LF}$ and ()$^{RF}$ refer to left fold point and right fold point, respectively. Indeed, in the $(a, b, c, \delta)$-space, the two points $p_1^{LF} = (a_1^{LF}, 0, c_1^{LF}, \delta^{LF})$ and $p_1^{RF} = (a_1^{RF}, 0, c_1^{RF}, \delta_1^{RF})$ (where $a_1^{LF}, a_1^{RF}, \delta_1^{LF}$ and $\delta_1^{RF}$ are obtained from $c_1^{LF}$ and $c_1^{RF}$ using Eqs. (C.41a) and (C.42a) and the b coordinate is null by Eq. (C.38a)) are generally called fold points.

Similarly, the local extrema of $H_2(c)$, denoted as $c_2^{LF}$ and $c_2^{RF}$, are obtained (again the explicit expressions are given by Eq. (C.67) in Appendix C.8) and exist if the following relations hold

$$ \beta > \frac{1}{2} \qquad (C.45) $$

and

$$ \mu < \frac{2\beta - 1}{\sqrt{3}}. \qquad (C.46) $$
<figure><img src="image_88.png" /><figcaption>**FIGURE C.6.** Typical example of the critical manifold. (a) In the (c, a)-plane. The branches $a = H_1(c)$ ($M_{01}$) and $a = H_2(c)$ ($M_{02}$) are plotted in red and green, respectively. The fold points $(c_i^{LF}, a_i^{LF})$ and $(c_i^{RF}, a_i^{RF})$ ($i = 1, 2$) (blue points), the points $(c_i^U, a_i^{LF})$ and $(c_i^D, a_i^{RF})$ (gray points) and the point $(c^I, a^I)$ (with $a^I = H_1(c^I) = H_2(c^I)$, see Eq. (C.39)) (magenta point) are also represented. A zoom is made to better see the shape of $M_{02}$. (b) In the (b, c, a)-space. (c) Zoom of (b) corresponding to the orange colored rectangular cuboid area. The following parameters are used: $\mu = 0.1, \beta = 0.65$ and $\alpha = 2$.</figcaption></figure>

If the conditions (C.45) and (C.46) are not satisfied, the function $H_2$ no longer has local extrema. As above, for the branch $M_{01}$, one can define left and right fold points for the branch $M_{02}$ as $p_2^{LF} = (a_2^{LF}, b_2^{LF}, c_2^{LF}, \delta^{LF})$ and $p_2^{RF} = (a_2^{RF}, b_2^{RF}, c_2^{RF}, \delta_2^{RF})$ (where $a_2^{LF}$, $a_2^{RF}$, $b_2^{LF}$, $b_2^{RF}$, $\delta_2^{LF}$ and $\delta_2^{RF}$ are obtained from $c_2^{LF}$ and $c_2^{RF}$ using Eqs. (C.38b), (C.41b) and (C.42b)).

Other important scalars for the present study, denoted as $c_1^D$ and $c_1^U$, are obtained solving $H_1(c_1^{RF}) = H_1(c_1^D)$ and $H_1(c_1^{LF}) = H_1(c_1^U)$. Similarly, if $c_2^{RF}$ and $c_2^{LF}$ exist, solving $H_2(c_2^{RF}) = H_2(c_2^D)$ and $H_2(c_2^{LF}) = H_2(c_2^U)$ yields $c_2^D$ and $c_2^U$. The expressions of $c_1^D, c_1^U, c_2^D$ and $c_2^U$ are given by Eqs. (C.68) and (C.69).

A typical example of the critical manifold in the (c, a)-plane is shown in Fig. C.6(a). The fold points $(c_i^{LF}, a_i^{LF})$ and $(c_i^{RF}, a_i^{RF})$ ($i = 1, 2$), the points $(c_i^U, a_i^{LF})$ and $(c_i^D, a_i^{RF})$ and the point $(c^I, a^I)$ (with $a^I = H_1(c^I) = H_2(c^I)$, see Eq. (C.39)) are represented. Moreover, a zoom has been made to better see the shape of $M_{02}$. The same example of the critical manifold in the (b, c, a)-space is shown in Fig. C.6(b).

Because in Fig. C.6 the conditions (C.44), (C.45) and (C.46) are satisfied, we can see that both $M_{01}$ and $M_{02}$ are S-shaped which suggests that several scenarios of relaxation oscillations are possible. To describe them (and other possible regimes), that requires knowing the stability of the critical manifold as well as the fixed points of the APMD. This is done in the next sections.

### C.4.2 Stability analysis of the critical manifold

The stability of the critical manifold $M_0$ is now investigated. For that, the Jacobian matrix $J_g$ of Eqs. (C.31b) to (C.31d) is computed and evaluated at any point of each branch of the critical manifold.
We obtain

$$ \mathbf{J}_g(p_1) = \begin{pmatrix} -H_{01}(c)^2 & 0 & 0 \\ \frac{\partial H_{02}}{\partial b}(0, c) & \frac{\partial H_{02}}{\partial c}(0, c) & -cH_{03}(0, c) \\ \frac{\partial H_{03}}{\partial b}(0, c) & \frac{\partial H_{03}}{\partial c}(0, c) + \frac{H_{03}(0, c)}{c} & \frac{H_{02}(0, c)}{c} \end{pmatrix}, \qquad (C.47) $$

for a point $p_1 \in \mathcal{M}_{01}$, and

$$ \mathbf{J}_g(p_2) = \begin{pmatrix} 2H_{01}(c)^2 & -2H_{01}(c)^2 H'_{01}(c) & 0 \\ \frac{\partial H_{02}}{\partial b}(H_{01}(c), c) & \frac{\partial H_{02}}{\partial c}(H_{01}(c), c) & -cH_{03}(H_{01}(c), c) \\ \frac{\partial H_{03}}{\partial b}(H_{01}(c), c) & \frac{\partial H_{03}}{\partial c}(H_{01}(c), c) + \frac{H_{03}(H_{01}(c), c)}{c} & \frac{H_{02}(H_{01}(c), c)}{c} \end{pmatrix}, \qquad (C.48) $$

for a point $p_2 \in \mathcal{M}_{02}$.

The critical manifold $\mathcal{M}_0$ contains both regular points defined as $\det(\mathbf{J}_g(p_i)) \neq 0$ ($i = 1, 2$) and singular points characterized by $\det(\mathbf{J}_g(p_i)) = 0$. Depending on the eigenvalues of $\mathbf{J}_g(p_i)$, the set of regular points can be split into several parts (see definitions in Chap. 3 of [40]), including normally hyperbolic *attracting* parts (all the eigenvalues have strictly negative real parts), normally hyperbolic *repelling* parts (all the eigenvalues have strictly positive real parts), normally hyperbolic saddle type parts (neither attracting nor repelling) and non normally hyperbolic parts (at least one eigenvalue has a zero real part). Note that at singular points the critical manifold is also not normally hyperbolic. In previous definitions, *normally* means that each point of $\mathcal{M}_0$ must be hyperbolic only in the the direction normal (i.e., non tangent) to itself (see e.g. Definition 2.3.4 and the explanations below [40]).

It can be shown that

$$ \det(\mathbf{J}_g(p_1)) = -H_{01}(c)^2 \frac{dH_1}{dc}(c) \frac{\sqrt{k_1^2 + k_2^2} \sqrt{c^2 H_{03}(0, c)^2 + H_{02}(0, c)^2}}{c} \qquad (C.49) $$

and

$$ \det(\mathbf{J}_g(p_2)) = H_{01}(c)^2 \frac{dH_2}{dc}(c) \frac{2\sqrt{k_1^2 + k_2^2} \sqrt{c^2 H_{03}(H_{01}(c), c)^2 + H_{02}(H_{01}(c), c)^2}}{c}. \qquad (C.50) $$

From Eqs. (C.49) and (C.50) we can deduce that the singular points of the critical manifold are, if they exist, the left and right fold points ($c_i^{\text{LF}}, a_i^{\text{LF}}$) and ($c_i^{\text{RF}}, a_i^{\text{RF}}$) ($i = 1, 2$) of each branch $\mathcal{M}_{0i}$ (the blue points in Fig. C.6) and also the point of intersection ($c^i, a^i$) between these two branches (the magenta point in Fig. C.6).

The three eigenvalues of $\mathbf{J}_g(p_1)$ (resp. $\mathbf{J}_g(p_2)$) are numerically computed; real and imaginary parts of these eigenvalues are plotted in Fig. C.7(a) (resp. Fig. C.7(b)). Comparing Figs. C.6 and C.7, we can see that the fold points and the point of intersection between the two branches of the critical manifold are actually points for which the latter is not normally hyperbolic. In addition, we can observe that the branch $\mathcal{M}_{02}$ loses its normal hyperbolicity at a regular point for $c \approx 0.33$ (see Fig. C.7(b) top). This corresponds to a Hopf bifurcation since at this point an equilibrium of the fast subsystem loses stability, as a pair of complex conjugate eigenvalues crosses the complex plane imaginary axis. Using the Routh-Hurwitz stability criterion and the Cardano's method, the analytical expression of this regular non normally hyperbolic fixed point of the branch $\mathcal{M}_{02}$ of the critical manifold can be obtained. The corresponding expression of $c$, denoted by $c^{\text{FP}}$, is

$$ c^{\text{FP}} = \sqrt{x_k} = \sqrt{2\sqrt{-\frac{p}{3}}\cos\left(\frac{1}{3}\arccos\left(\frac{3q}{2p}\sqrt{\frac{3}{-p}}\right) + \frac{2\pi}{3}\right) - \frac{b_2}{3b_1}} \qquad (C.51) $$

The definitions of $p$ and $q$ and details on the computation of $c^{\text{FP}}$ are given in Appendix C.9. As previously, the other coordinates of the fixed point, i.e., $b^{\text{FP}}$, $a^{\text{FP}}$ and $\delta^{\text{FP}}$, can be deduced from $c^{\text{FP}}$ using Eqs. (C.38b), (C.41b) and (C.42b)), respectively.

The previous stability analysis is summarized in Fig. C.8 in which each branch $\mathcal{M}_{01}$ and $\mathcal{M}_{02}$ of the critical manifold is decomposed as follows

$$ \mathcal{M}_{01} = \mathcal{M}_{01}^{\text{st},1} \cup p^{\text{I}} \cup \mathcal{M}_{01}^{\text{st},2} \cup p_1^{\text{LF}} \cup \mathcal{M}_{01}^{\text{st},2} \cup p_1^{\text{RF}} \cup \mathcal{M}_{01}^{\text{a},2} \qquad (C.52) $$

with $p^I = (a^I, 0, c^I, \delta^I)$ (where $a^I$ and $\delta^I$ are obtained from $c^I$ using Eqs. (C.41a) and (C.42a) and the $b$ coordinate is null by Eq. (C.38a)), and

$$ \mathcal{M}_{02} = \mathcal{M}_{02}^{\text{a},1} \cup p_2^{\text{LF}} \cup \mathcal{M}_{02}^{\text{st},1} \cup p_1^{\text{RF}} \cup \mathcal{M}_{01}^{\text{a},2} \cup p^{\text{FP}} \cup \mathcal{M}_{02}^{\text{st},2} \cup p^{\text{I}} \qquad (C.53) $$
<figure><img src="image_89.png" /><figcaption>**FIGURE C.7.** The three eigenvalues of (a) $J_g(p_1)$ given by (C.47) and (b) $J_g(p_2)$ given by (C.48). Top: real parts; Bottom: imaginary parts. The eigenvalues are numerically computed using the same BNES parameters as in Fig. C.6.</figcaption></figure>

with $p^{\text{FP}} = (a^{\text{FP}}, b^{\text{FP}}, c^{\text{FP}}, \delta^{\text{FP}})$ (where $a^{\text{FP}}, b^{\text{FP}}$ and $\delta^{\text{FP}}$ are obtained from $c^{\text{FP}}$ using Eqs. (C.38b), (C.41b) and (C.42b)). In both previous expressions the superscripts “st” and “a” indicate the saddle-type and attracting natures of the critical manifold, respectively.

The structure of $M_{02}$ can be modified if (C.45) or (C.46) does not hold, in this case $M_{02}$ is not S-shaped anymore and $M_{02}^{\text{st},1}$ disappears. Moreover, if $p^{\text{FP}}$ joins $p^I$, the branch $M_{02}^{\text{st},2}$ disappears.

From the previous analysis, we can already give a partial description of the APMD behavior. In phase space, if the trajectory is not in the vicinity of the critical manifold, it evolves fast to an attracting part of the critical manifold. These fast epochs are approximately described by the fast subsystem (C.31). Bistability has been observed since several points of the critical manifold can be stable for a given set of parameters and a given value of the slow variable $a$ (which, due to (C.31a), is a parameter for the fast subsystem). In the present paper the basins of attraction of each attracting branch of the critical manifold are not rigorously investigated, together with the possible other kinds of solutions of the fast subsystem (such as periodic, quasiperiodic or even chaotic motions). This can be investigated in the future.

In the vicinity of the critical manifold, the trajectory evolves slowly. These slow epochs are approximately described by the slow subsystem (C.30), its analysis is performed in the next section.
<figure><img src="image_90.png" /><figcaption>FIGURE C.8. The different parts of the critical manifold $M_0$. (a) In the (c, a)-plane. For the branch $M_{01}$ (in red) one has (see Eq. (C.52)) from left to right: the first saddle-type branch $M_{01}^{st,1}$, the first attracting branch $M_{01}^{a,1}$, the second saddle-type branch $M_{01}^{st,2}$ and the second attracting branch $M_{01}^{a,2}$. These branches are respectively connected to each other by the point $p^I$ (magenta point) and the left and right fold points $p_1^{LF}$ and $p_1^{RF}$ (blue points). For the branch $M_{02}$ (in green) one has (see Eq. (C.53)) from left to right: the first attracting branch $M_{02}^{a,1}$, the first saddle-type branch $M_{02}^{st,1}$, the second attracting branch $M_{02}^{a,2}$ and the second saddle-type branch $M_{02}^{a,2}$ which ends at the point $p^I$. These branches are respectively connected to each other by the left and right fold points $p_2^{LF}$ and $p_2^{RF}$ (blue points) and the point $p^{FP}$ (green point). (b) In the (b, c, a)-space. (c) Zoom of (b) corresponding to the orange colored rectangular cuboid area.</figcaption></figure>

### C.4.3 Asymptotic analytical expression of the fixed points of the amplitude-phase modulation dynamics

As mentioned previously, the fixed points of the APMD (C.29) (or (C.27)) characterize periodic solutions of the original system (E.3). Since $0 < \epsilon \ll 1$, these fixed points can be asymptotically approximated by those of the slow subsystem (C.30). Moreover, to obtain a simple analytical expression of these fixed points, it is first necessary to linearize the function $f$ (see Eq. (C.28a)) around $a = 0$.

Then, using successively Eqs. (C.38a), (C.41a) and (C.42a) and Eqs. (C.38b), (C.41b) and (C.42b) the slow dynamics (C.30) on the branches $M_{01}$ and $M_{02}$ can be considered only with respect to $c$ as

$$ \frac{dH_1}{dc}(c)c' = F_1(c) \qquad (C.54) $$

and

$$ \frac{dH_2}{dc}(c)c' = F_2(c), \qquad (C.55) $$

respectively. The expressions of the functions $F_1$ and $F_2$ in the previous equations are given in Appendix C.8.

Consequently, the desired fixed points are solutions of $F_1(c) = 0$ and $F_2(c) = 0$ of which $c = 0$ is a trivial solution. In the (b, c, a)-space, on $M_{01}$ this gives (0, 0, 0) which correspond to $p_0^{st}$ (see Eq. (C.3)). On $M_{02}$ that yields $(0, \pm\sqrt{\beta/\alpha}, 0)$ (see Eq. (C.38b)) which correspond to $p_1^{st}$ and $p_2^{st}$ (see Eq. (C.4)).

Then the non trivial solutions of $F_1(c) = 0$ and $F_2(c) = 0$ are computed. The terms in parentheses in the numerator of $F_1$ and $F_2$ are quadratic polynomials with respect to $c^2$. The corresponding positive solutions in
terms of the variable c are $c_{1,1}^*$ and $c_{1,2}^*$ obtained solving $F_1(c) = 0$, and $c_{2,1}^*$ and $c_{2,2}^*$ obtained solving $F_2(c) = 0$. The explicit expressions of $c_{1,1}^*$, $c_{1,2}^*$, $c_{2,1}^*$ and $c_{2,2}^*$ are given by Eqs. (C.71) and (C.72) in Appendix C.8. One can easily verify from Eqs. (C.71) and (C.72) that $c_{2,1}^* < c_{2,2}^* < c_{1,1}^* < c_{1,2}^*$ and that these fixed points no longer exist when they become complex for

$$ \sigma > \sigma^S = \frac{1}{\mu}. \qquad (C.56) $$

Other sets of parameters can make the above expressions complex, except $c_{1,2}^*$ which is always real if Eq. (D.34) holds.

Finally, the stability of the fixed points $c_{1,i}^*$ and $c_{2,i}^*$ ($i = 1, 2$) is obtained checking the sign of $\frac{d}{dc} \left( \frac{F_1(c)}{\frac{dH_1}{dc}(c)} \right) \bigg|_{c=c_{1,i}^*}$ and $\frac{d}{dc} \left( \frac{F_2(c)}{\frac{dH_2}{dc}(c)} \right) \bigg|_{c=c_{2,i}^*}$ ($i = 1, 2$), respectively.

The analysis presented in this section provides an approximate description of the slow dynamics of the APMD in the vicinity of an attracting branch of the critical manifold. However, slow dynamics around families of periodic, quasiperiodic or chaotic solutions of the fast subsystem is not investigated. Studies inspired from the work of Berglund and Gentz [22] on invariant manifold tracking the family of periodic orbits could be again the subject of future works.

## C.5 ASYMPTOTIC GLOBAL STABILITY ANALYSIS OF THE AMPLITUDE-PHASE MODULATION DYNAMICS

Knowing the fixed points of the APMD and their stability, obtained in Section C.4.3, together with the analysis of the critical manifold presented in Section ??, allow to perform the global stability analysis of the APMD, again asymptotically in the case of $\epsilon \to 0$. This will enable us to interpret the observations made in Fig. C.4. The analysis is possible analytically because the critical manifold is one-dimensional. Therefore the slow dynamics is also one-dimensional and described by Eqs. (C.54) and (C.55).

From Eqs. (C.11), (C.16), (C.17), (C.20), (C.21) and the rules $b = b_0 + \epsilon b_1$ and $C = C_0 + \epsilon C_1$, one has

$$ y(t) = - \frac{b(t)\epsilon + b(t) - a(t)\cos(\theta(t) + t) + c(t)\cos(t + \varphi(t)) + \frac{1}{32}a(t)^3\lambda\epsilon\sin(3(\theta(t) + t))}{\epsilon + 1}. \qquad (C.57) $$

For comparison purposes, the time series $y(t)$ reconstructed from numerical simulations of the APMD using Eq. (C.57) are shown in Fig. C.9. In Figs. C.9(a) to C.9(e) the parameters used are the same as in Figs. C.4(a) to Figs. C.4(e). We can see that the APMD can reproduce regimes of the same nature as those obtained with the full order system (E.3). In Figs. C.9(f) and C.9(g), even if the APMD can reproduce these two chaotic regimes, they are here obtained for smaller values of the bifurcation parameter of $\sigma$, namely $\sigma = 1.22$ and 0.77 for Figs. C.9(f) and C.9(g), respectively. In general when chaotic regimes occur the mitigation limit observed on the numerical simulation of the APMD is smaller than that obtained with the full order system (E.3).

To explain the nature of the different regimes shown in Figs. C.4 and C.9, the trajectories of the APMD in the (c, a)-plane (obtained from the same simulations as in Fig. C.9) are superimposed on the critical manifold given by Eq. (C.41) in Fig. C.10. The structure of the critical manifold (i.e, stability, fold points ($c_i^{LF}$, $a_i^{LF}$) and ($c_i^{RF}$, $a_i^{RF}$) ($i = 1, 2$), points ($c_i^U$, $a_i^{LF}$) and ($c_i^D$, $a_i^{RF}$) and the point ($c^I$, $a^I$)) is depicted in the same way as in Figs. C.6 and C.8. Fig. C.10 shows also the fixed points of the APMD obtained in Section C.4.3, i.e. the trivial fixed points 4 and the non trivial fixed points ($c_{1,i}^*$, $H_1(c_{1,i}^*)$) (in black color) and ($c_{2,i}^*$, $H_2(c_{2,i}^*)$) (in orange color) with $i = 1, 2$ and $c_{1,i}^*$ and $c_{2,i}^*$ given by Eqs. (C.72) and (C.72), respectively. These fixed points are depicted by a ⋆ when they are stable and a ● when they are unstable.

The different regimes can now be explained as follows :

* Fig. C.10(a). The trivial fixed point on $M_{02}^{a,1}$ is stable (depicted by an orange ⋆). The other fixed points on $M_{02}^{a,1}$ is unstable (depicted by an orange ●). As long as the initial condition causes the APMD trajectory to arrive (after a fast epoch) on $M_{02}^{a,1}$ between the unstable fixed point and the stable fixed point, the latter is reached after a slow epoch along $M_{02}^{a,1}$, as in Fig. C.10(a). The initial condition in the zoom (orange frame in Fig. C.10(a)) is depicted by a blue ●.

* Fig. C.10(b). There is only the unstable trivial fixed point on $M_{02}^{a,1}$ (for this set of parameters one has $M_{02} = M_{02}^{a,1}$). The trajectory follows slowly $M_{02}^{a,1}$ until it disappears at ($c^I$, $a^I$) and then follows $M_{01}^{a,1}$.

<footnote n="4">4. Remember that they correspond to $p_0^{st}$ when they are on $M_{01}$ and to $p_1^{st}$ (or $p_2^{st}$) when they are on $M_{02}$. They are superimposed in the (c, a)-plane but not in the (b, c, a)-space (see Section C.4.3).</footnote>
<figure><img src="image_91.png" /><figcaption>**FIGURE C.9.** Time series $y(t)$ reconstructed from numerical simulations of the APMD using Eq. (C.57). Equilibrium positions $y_1^{st}$ (dashed green line) and $y_2^{st}$ (dashed red line). From (a) to (e) the parameters used are the same as in Figs. C.4(a) to Figs. C.4(e). In (f) and (g) the value of $\sigma$ is changed compared to Figs. C.4(f) and C.4(g), namely $\sigma = 1.22$ and $0.77$, respectively.</figcaption></figure>
slowly again. Finally, the trajectory stops on a stable fixed point (depicted by a black ⋆). Because on $M_{01}$ $b=0$ (see Fig. C.8(b)), this scenario of the APMD explains the periodic regime PR1 observed in Fig. C.4(b).

* Fig. C.10(c). The fixed point on $M_{01}^{a,1}$ moves to the saddle-type part $M_{01}^{st,2}$. Consequently, the trajectory evolves slowly to the left fold point $(c_1^{LF}, a_1^{LF})$ and undergoes a fast jump to $M_{01}^{a,2}$. Because the unstable point on $M_{01}^{a,2}$ (the black ●) is above $(c_i^U, a_i^{LF})$ (the gray ■) we then observe a succession of slow evolutions near the attracting branch of $M_{01}$ and fast jumps (from $(c_1^{LF}, a_1^{LF})$ to $(c_i^U, a_i^{LF})$ and from $(c_1^{RF}, a_1^{RF})$ to $(c_i^D, a_i^{RF})$). These so-called *relaxation oscillations* on $M_{01}$ explain the SMR1 observed in Fig. C.4(c).

* Fig. C.10(d). There is a stable non trivial fixed point on $M_{02}^{a,1}$ (the orange ⋆). The trivial fixed point and the other non trivial fixed point on $M_{02}^{a,1}$ are unstable (the two orange ●). As long as the initial condition causes the APMD trajectory to arrive on $M_{02}^{a,1}$ below the unstable non trivial fixed point, the stable fixed point is reached after a slow epoch along $M_{02}^{a,1}$, as in Fig. C.10(c). This scenario for the APMD corresponds to the periodic intra-well motion PR1 observed in Fig. C.4(d).

* Fig. C.10(e). A scenario of relaxation oscillations on $M_{02}$ (similar to that observed on $M_{01}$ in Fig. C.10(c)) is observed in agreement with the intra-well SMR2 of Fig. C.4(d). Again, as long as the initial condition causes the APMD trajectory to arrive on $M_{02}^{a,2}$ below the unstable non trivial fixed point, these relaxation oscillations will occur.

* Fig. C.10(f). For the sake of readability of the figure, only steady state motion is shown. This scenario is also illustrated in Fig. C.11(a), showing the time series of the variables $b(t)$, $c(t)$ and $y(t)$ (computed again using Eq. (C.57)), and in Fig. C.12(a) depicting the trajectory of the APMD in the $(b, c, a)$-space. For this scenario and the next one, when we speak of an equilibrium (resp. a periodic) regime during a slow epoch of the APMD, this means that for the set of parameters and the value of $a$ considered, the fast subsystem has a stable equilibrium (resp. periodic) solution. The APMD therefore follows an invariant manifold tracking these equilibrium (resp. periodic) solutions and parameterized by the slow variable $a$. For the full order dynamics (referred to as FOD in Fig. C.11(a)), which includes oscillations at the natural frequency of the VdP oscillator, that corresponds to invariant manifold tracking of periodic (resp. quasiperiodic solutions). Chaotic motions are also encountered. The successive phases of the motion are described in Fig. C.11(a). During a complete cycle, the trajectory of the APMD :

- Follows $M_{02}^{a,1}$ (SE1: APMD, Equilibrium 1; FOD, Periodic 1). For illustration purposes, results of the numerical simulation of the fast subsystem (C.31) for $a = 0.07$ are shown in Fig. C.13(a). We can see that with suitably chosen initial conditions, the fast subsystem dynamics reaches an equilibrium which corresponds to a point of $M_{02}^{a,1}$.

- Undergoes a fast jump (FE2) from $(c_2^{LF}, a_2^{LF})$ and the corresponding left branch of $M_{02}$ is saddle-type (i.e., $M_{02}^{st,2}$).

- Is periodic (Periodic 1) and quasiperiodic (Quasiperiodic 1) for the FOD precisely because $M_{02}$ is saddle-type and then chaotic (Chaotic 1). For the APMD (resp. FOD), the phases Periodic 1 (resp. Quasiperiodic 1) and Chaotic 1 constitute the SE2. Numerical solutions of the fast subsystem (C.31) for $a = 0.14$ and $0.25$ are shown in Figs. C.13(b) and C.13(c), respectively. Again with suitable initial conditions, the fast subsystem dynamics reaches a periodic solution when $a = 0.14$ and a chaotic attractor when $a = 0.25$, in agreement with what we observe for the APMD during Periodic 1 and Chaotic 1 phases, respectively (see the observed APMD for the same values of $a$ in Fig. C.12(a)).

- Undergoes a fast jump (FE3) to $M_{01}^{a,2}$.

- Follows $M_{01}^{a,2}$ (SE3: APMD, Equilibrium 2; FOD, Periodic 2). Results of the numerical simulation of the fast subsystem (C.31) are shown in Fig. C.13(a) for $a = 0.4$. We can see here that, after any transient dynamics has died out, the fast subsystem dynamics reaches an equilibrium which corresponds to a point of $M_{01}^{a,2}$.

- Finally undergoes a fast jump (FE1') from $(c_1^{RF}, a_1^{RF})$ to $M_{02}^{a,1}$, and so on.

At each cycle there is the possibility of reaching one or other of the opposite $M_{02}$ branches (see Fig. C.12(a)). Concerning the start of the chaotic phases, we do not know whether this is really a chaotic movement of the fast subsystem or a quasi-periodic solution. Again, this may be subject of future work with a more thorough bifurcation analysis of the fast subsystem, including the determining of basins of attraction in case of multistability. This scenario of the APMD is in agreement with what
we observed in Fig. C.5(a) except that the Periodic 1 and Quasiperiodic 1 phases (for the FOD) appear to be reversed.

* Fig. C.10(g). Again, only steady state motion is shown and the scenario is also illustrated in Fig. C.11(b), showing the time series of the variables *b*(*t*), *c*(*t*) and *y*(*t*) (computed again using Eq. (C.57)), and in Fig. C.12(b) depicting the trajectory of the APMD in the (*b*, *c*, *a*)-space. The successive phases of the motion are described in Fig. C.11(b). When we observe the APMD, we see a periodic and then a chaotic regime during phases 1 and 1'. It is difficult to differentiate between these two phases on the *y*(*t*) signal from the FOD. However, the latter is qualitatively equivalent to what we observed in Fig. C.5(b). Numerical solutions of the fast subsystem (C.31) for *a* = 0.22, 0.35 and 0.45 are shown in Figs. C.14(a), C.14(b) and C.14(c), respectively. In agreement with what we observe for the APMD during Periodic 1, Chaotic 1 and Equilibrium phases (see the observed APMD for the same values of *a* in Fig. C.12(b)), the fast subsystem dynamics reaches a periodic regime, a chaotic motion and an equilibrium for *a* = 0.22, 0.35 and 0.45, respectively.

In each of the previous situations we can see that if there is a passage on $M_{01}$, it is the largest unstable fixed point (black • on $M_{01}^{a,2}$) which prevents the system from tipping over into a harmful situation. If the trajectory ends up on $M_{01}^{a,2}$, above this fixed point, it runs along $M_{01}^{a,2}$ towards infinity (i.e. harmful situation) when the linearized VdP is considered (see Section C.4.3). If the nonlinear VdP is considered, the trajectory reaches a stable fixed point with large amplitude similar to that of the VdP alone (i.e. also harmful situation).

Of course, the global stability analysis presented here is only partial. In the future, it will be necessary to rigorously research the possible solutions of the fast subsystem, as well as their stability. It will then be necessary to find the invariant manifolds of the APMD tracking these solutions. However, this study enables us to interpret a certain number of regimes observed on numerical simulations.

## C.6 CONCLUSION

In this paper a Van der Pol oscillator coupled to a bistable nonlinear energy sink BNES) has been studied. Numerical simulations have shown that this system can undergo a multitude of motions including different types of periodic regimes and so-called strongly modulated responses (SMR) as well as chaotic regimes. We also show that a BNES can be much more efficient than a classical cubic NES but this is not robust since a little perturbation can switch the system from harmless to harmful situations. However, even in the most unfavorable cases, we were able to find a set of parameters for which the BNES performs better than the NES.

Then, in an original way, the so-called Multiple Scale/Harmonic Balance Method (MSHBM) is used and modified (compared to its usual use) to consider the specific feature of the BNES, i.e., that it can have a nonzero-mean oscillating motion. This allowed us to derive an amplitude-phase modulation dynamics (APMD) which can reproduce the complex behavior of the initial system.

Because of the presence of a small perturbation parameter (i.e., the mass ratio between the BNES and the VdP oscillator), the APMD is governed by two different time scales. More precisely, in its real form it appears as a (3,1)-fast-slow system. The motion of such a fast-slow system consists in a succession of slow and fast epochs. In the case of a classic cubic NES, the APMD is a (2,1)-fast-slow system whose behavior is now well understood. For example, the most complex responses are relaxation oscillations leading to so-called strongly modulated responses (SMR) for the initial system. Here, due to its (3,1)-fast-slow nature, the APMD can undergo a multitude of different motions including periodic and quasiperiodic regimes of different types and even chaotic regimes, in agreement with what has been observed on numerical simulations of the initial system. The fast-slow analysis of the APMD has been conducted within the framework of the geometric singular perturbation theory. By the computation of the so-called critical manifold and the analytical expressions of the APMD fixed points, a global stability analysis has been performed. This enabled us to interpret a certain number of regimes observed on numerical simulations of the initial system.

The methodology proposed here is applied to the study of passive mitigation of self-sustained oscillations using a BNES. In this context, the Van der Pol oscillator is used here as an archetypal self-sustained oscillator. However, in mechanical systems which are more than one-dimensional, self-sustained oscillations appear in general because of mode-coupling instabilities. This is the case, for example, with the flutter instability of an aircraft wing or with friction-induced vibrations in braking systems. Classical cubic NES have demonstrated their ability to mitigate this type of self-sustained oscillations (see e.g. [a8, a14, 78]). The present results suggest that BNES could be even more effective. Moreover, if the primary mechanical structure has only one unstable mode (i.e., the modal analysis reveals that only one mode can have an eigenvalue with a positive real part) the methodology proposed in [a8, 78] shows how to reduce the primary structure dynamics on its unstable mode. This allows to analyze the system as in the case of a one-dimensional primary structure. These
<figure><img src="image_92.png" /><figcaption>**FIGURE C.10.** Trajectories of the APMD in the (c, a)-plane (obtained from the same simulations as in Fig. C.9) are superimposed on the critical manifold given by Eq. (C.41). The structure of the critical manifold (i.e, stability, fold points ($c_i^{LF}$, $a_i^{LF}$) and ($c_i^{RF}$, $a_i^{RF}$) ($i = 1, 2$), points ($c_i^U$, $a_i^{LF}$) and ($c_i^D$, $a_i^{RF}$) and the point ($c^I$, $a^I$)) is depicted in the same way as in Figs. C.6 and C.8. Fixed points of the APMD obtained in Section C.4.3, i.e. the trivial fixed points and the non trivial fixed points ($c_{1,i}^*$, $H_1(c_{1,i}^*)$) (in black color) and ($c_{2,i}^*$, $H_2(c_{2,i}^*)$) (in orange color) with $i = 1, 2$ and $c_{1,i}^*$ and $c_{2,i}^*$ given by Eqs. (C.71) and (C.72), respectively. These fixed points are depicted by a ⋆ when they are stable and a ● when they are unstable.</figcaption></figure>
<figure><img src="image_93.png" /><figcaption>**FIGURE C.11.** Same numerical simulations as in (a) Figs. C.10(f) and (b) C.10(g). For both (a) and (b) from top to bottom one has the time series *b*(*t*) (gray), (*c*) (gray) and *y*(*t*) (blue). The equilibrium positions *y*<sub>1</sub><sup>st</sup> (dashed green line) and *y*<sub>2</sub><sup>st</sup> (dashed red line) are depicted in the plot of *b*(*t*). One has APMD for amplitude-phase modulation dynamics, FOD for full order dynamics, FE for fast epoch and SE for slow epoch.</figcaption></figure>