* Err is now the relative error of the optimal design value $\varphi_h^{\text{opt}}$ obtained with the ME-gPC method, compared to that obtained with the reference method :

$$Err = \frac{|(\varphi_h^{\text{opt}})_{\text{ME-gPC}} - (\varphi_h^{\text{opt}})_{\text{ref}}|}{(\varphi_h^{\text{mid}})_{\text{ref}}}. \qquad (\text{F.24})$$

### F.6.2.1 Reference optimization

The same strategy has been applied for the reference optimization of the NES cubic stiffness $\varphi_h$. First, to get an idea of the shape of the PHSSR as a function of $\varphi_h$, a small number of simulations is performed with a discretization step equal to 10 on the design space $\varphi_h \in [0, 900 000] \text{ N}\cdot\text{kg}^{-1}\cdot\text{m}^{-3}$. For each value of $\varphi_h$ the PHSSR is computed using 1000 values of $\mu$ within the interval $[0, 0.4]$. Fig. F.9 shows the resulting evolution of the PHSSR as a function of $\varphi_h$. The mitigation limit appears to decrease monotonically with the cubic stiffness value. The optimal value $\varphi_h^{\text{opt}}$ is therefore the lowest value of $\varphi_h$ for which the jump criterion (F.16) is met and is found equal to 100 000 N\cdotkg^{-1}\cdotm^{-3}.

The next step then consists in focusing on a reduced design space composed of the lowest values of $\varphi_h$ in order to understand the dynamic behavior of the system for these values and to refine the optimization result. The retained design interval $\varphi_h \in [0, 150 000] \text{ N}\cdot\text{kg}^{-1}\cdot\text{m}^{-3}$ is uniformly discretized into 100 values, and for each one the mitigation limit is directly computed from (F.2) using again 1000 values of $\mu$ within $[0, 0.4]$ interval. The number of simulations is thus equal to $100 \times 1000 = 100 000$.

<figure><img src="image_123.png" /><figcaption>**FIGURE F.9.** PHSSR as a function of $\varphi_h$ for $\eta_h = 0.024$ and $\varphi_h \in [0, 900 000] \text{ N}\cdot\text{kg}^{-1}\cdot\text{m}^{-3}$. The mitigation limits are assessed using the reference method.</figcaption></figure>

Fig. F.10 shows the evolutions of the jump criterion (F.16) and of the PHSSR as functions of $\varphi_h$, respectively. For low values of $\varphi_h$ ($\varphi_h \in [0, 46 970] \text{ N}\cdot\text{kg}^{-1}\cdot\text{m}^{-3}$ a jump exists but, as mentioned before, we consider that the difference between the amplitudes of the limit cycles before and after the jump is too small for the concerned regimes to be qualified as mitigated. The difficulty of this optimization problem is then to find the lowest value of $\varphi_h$ for which the discontinuity jump is high enough. It is clear that the solution depends on the value of the threshold for the jump criterion (here 50%). The optimal value here is therefore $(\varphi_h^{\text{opt}})_{\text{ref}} = 46 970 \pm 3.22\% \text{ N}\cdot\text{kg}^{-1}\cdot\text{m}^{-3}$, which corresponds to the maximum PHSSR value $\text{PHSSR}^{\text{max}} = 52.25\% \pm 0.19\%$. The PHSSR curve exhibits again stationary phases linked to the discretization step for $\mu$. Indeed, for two close values of $\varphi_h$, the discontinuities in the amplitude profiles as functions of $\mu$ appear for very close $\mu$ values whose difference is smaller than the discretization step.

### F.6.2.2 Optimization with the ME-gPC method

We denote by $(\varphi_h^{\text{opt}})_{\text{ME-gPC}}$ the optimal value of the NES cubic stiffness coefficient obtained by solving (F.17) where the amplitude $A(x_d, x_u)$ and the mitigation limit are assessed using the ME-gPC based method.
<figure><img src="image_124.png" /><figcaption>**FIGURE F.10.** Evolutions of the Jump Criterion (JC) and of the PHSSR as functions of $\phi_h$ obtained with the reference method.</figcaption></figure>

<table><caption>**TABLEAU F.2.** Comparison between the reference optimization and the ME-gPC based optimization considering deterministic design variables. The design parameter is $\phi_h$. The reference results are obtained as described in Sect. F.6.2.2.</caption><thead><tr><th></th><th>Nb of simulations</th><th>Computational time</th><th>PHSSR<sup>max</sup></th><th>E<sub>PHSSR<sup>max</sup></sub></th><th>φ<sub>h</sub><sup>opt</sup></th><th>Err</th></tr></thead><tbody><tr><td>Reference</td><td>100 000</td><td>1.49 days</td><td>52.25% ± 0.19%</td><td>—</td><td>46 970 ± 3.22%</td><td>—</td></tr><tr><td>ME-gPC<br>(θ<sub>1</sub> = 3.3 · 10<sup>-3</sup><br>and θ<sub>2</sub> = 1%)</td><td>6 652</td><td>2.64 hours</td><td>53.90% ± 1.86%</td><td>3.15%</td><td>33 333 ± 4.50%</td><td>29.03%</td></tr><tr><td>ME-gPC<br>(θ<sub>1</sub> = 4.2 · 10<sup>-3</sup><br>and θ<sub>2</sub> = 1%)</td><td>5 817</td><td>2.31 hours</td><td>53.13% ± 1.88%</td><td>1.68%</td><td>39 394 ± 3.40%</td><td>22.51%</td></tr><tr><td>ME-gPC<br>(θ<sub>1</sub> = 4.8 · 10<sup>-3</sup><br>and θ<sub>2</sub> = 1%)</td><td>5 706</td><td>2.25 hours</td><td>53.13% ± 1.88%</td><td>1.68%</td><td>45 454 ± 3.33%</td><td>3.22%</td></tr><tr><td>ME-gPC<br>(θ<sub>1</sub> = 3.5 · 10<sup>-3</sup><br>and<br>θ<sub>2</sub> = 0.1%)</td><td>7 023</td><td>2.96 hours</td><td>52.34% ± 0.19%</td><td>0.35%</td><td>46 970 ± 3.22%</td><td>0%</td></tr></tbody></table>

The gPC order is again set to $p = 1$ and the mitigation limit, as explained before, is approached by the upper bound of the last element found. Tab. F.2 compares the results obtained with this ME-gPC method, for several thresholds values $\theta_1$ and $\theta_2$, to those arising from the reference procedure (in which, as detailed in the previous section, the values of $\phi_h$ are taken from the design space [0, 150 000] N·kg<sup>-1</sup>·m<sup>-3</sup> which is uniformly discretized).

Fig. F.11 shows the evolution of the PHSSR as a function of $\phi_h$ for $\theta_1 = 3.5 \cdot 10^{-3}\%$ and $\theta_2 = 0.1\%$.

For a given value of the threshold $\theta_2$, it can be observed from Tab. F.2 that the precision of the results depends on the jump threshold $\theta_1$. Indeed, when the value of $\theta_1$ increases, the error PHSSR<sup>max</sup> decreases from 3.15% to 1.68% and the error Err decreases from 29.03% to 3.22%. Another benefit of increasing $\theta_1$ is that it reduces the number of required simulations, as can be seen in Tab. F.2. However, for higher values of $\theta_1$, the ME-gPC method appears unable to detect the discontinuity in the amplitude profile, as the jump in the amplitude values becomes lower than the threshold $\theta_1$, which makes the research of the optimum $\theta_1$ value difficult. Another way to reduce the error levels is to increase the value of the precision threshold $\theta_2$. The optimum value of the jump detection threshold $\theta_1$ depending on the precision threshold $\theta_2$, the optimum pair
of thresholds is finally found as $\theta_1 = 3.5 \cdot 10^{-3}$ and $\theta_2 = 0.1 \%$. The corresponding relative errors are 0.19 % for PHSSR$^{max}$ and $Err = 0 \%$ for $\varphi_h^{opt}$, while the number of simulations has increased from 5706 to 7023. It is therefore observed that the precision level depends on the accepted computational cost. In this study, the final computational cost (7 023) with the ME-gPC method remains much weaker than that of the reference method (100 000).

*Remark*: as in section F.6.1 for the optimization of the NES damping coefficient, identical values of PHSSR$^{max}$ may be observed with the ME-gPC method. The widths of the plateaus decrease when the precision threshold $\theta_1$ is reduced.

<figure><img src="image_125.png" /><figcaption>**FIGURE F.11.** PHSSR as a function of $\varphi_h$ for $\eta_h = 0.024$ and $\varphi_h \in [0 \ 150\ 000] \ \text{N} \cdot \text{kg}^{-1} \cdot \text{m}^{-3}$. The thresholds are respectively: $\theta_1 = 3.5 \cdot 10^{-3}$ and $\theta_2 = 0.1 \%$. The reference results are obtained as described in Sect. F.6.2.2.</figcaption></figure>

## F.7 RESULTS FOR STOCHASTIC DESIGN VARIABLES

The second optimization technique considers that the NES parameters are also uncertain with a known probability law. Thus, the design variables to be optimized are not directly the NES parameters but one of their statistics called hyperparameters. In this work, the design parameters follow uniform probability distributions and the considered hyperparameter is the mean $\tilde{x}_d$. Regarding the system parameters, as previously only the friction coefficient $\mu$ is assumed to be uncertain.

### F.7.1 Optimization of the NES damping coefficient

#### F.7.1.1 Reference optimization

In this section, the optimization problem (F.18) is solved with a single design variable (hyperparameter) defined as the mean $\hat{\eta}_h$ of the NES damping coefficient $\eta_h$, which is a stochastic parameter following a uniform law. The NES cubic stiffness is fixed and equal to $\varphi_h = 140\ 000 \ \text{N} \cdot \text{kg}^{-1} \cdot \text{m}^{-3}$. The design space is divided into $N_I = 5$ intervals as shown in Tab. F.3. For each interval the lower and upper bounds are determined as follows: $\eta_h^{\text{inf}} = \hat{\eta}_h - 0.1\hat{\eta}_h$ and $\eta_h^{\text{sup}} = \hat{\eta}_h + 0.1\hat{\eta}_h$ where $\hat{\eta}_h$ is the mean. The mean values are chosen so that the upper bound of one interval coincides with the lower bound of the next interval.

In each interval, $N_d = 10$ random samples of damping coefficient values $\eta_h$ are drawn according to the uniform probability law. For each of these 10 samples, the mitigation limit and the PHSSR are determined directly from the numerical simulation of (F.2) using 10000 values of the friction coefficient $\mu$. The objective function to be maximized is now the mean $E[\text{PHSSR}]$ of the 10 PHSSR values obtained in each interval.

Tab. F.3 shows the reference results of the $\hat{\eta}_h$ optimization. For each interval, the mean of the 10 samples of $\eta_h$ is calculated and compared to the theoretical $\hat{\eta}_h$ value in order to validate the number of samples. The mean PHSSR value $E[\text{PHSSR}]$ is also indicated, along with the error on $E[\text{PHSSR}]$ for a confidence level of 99%. The resulting values (errors lower than 0.15%) show that the number of samples $N_d = 10$ of $\eta_h$ in each interval is sufficient.
<table><caption>**TABLEAU F.3.** Results of the reference optimization (see Sect. F.7.1.1). The design parameter is $\hat{\eta}_h$.</caption><thead><tr><th>[η<sub>h</sub><sup>inf</sup>, η<sub>h</sub><sup>sup</sup>]</th><th>Theoretical mean η̂<sub>h</sub></th><th>Mean of the 10 samples of η<sub>h</sub></th><th>E[PHSSR] (%)</th><th>Error on E[PHSSR] for a confidence level = 99% (%)</th></tr></thead><tbody><tr><td>[0.0117,0.0143]</td><td>0.0130</td><td>0.0129</td><td>48.36</td><td>0.1417</td></tr><tr><td>[0.0143,0.0175]</td><td>0.0159</td><td>0.0161</td><td>49.31</td><td>0.1086</td></tr><tr><td>[0.0175,0.0214]</td><td>0.0194</td><td>0.0193</td><td>48.85</td><td>0.0220</td></tr><tr><td>[0.0214,0.0261]</td><td>0.0237</td><td>0.0234</td><td>50.07</td><td>0.0005</td></tr><tr><td>[0.0261,0.0319]</td><td>0.0290</td><td>0.0284</td><td>49.90</td><td>0.0277</td></tr></tbody></table>

<table><caption>**TABLEAU F.4.** Results of the optimization with the ME-gPC method. The design parameter is $\hat{\eta}_h$.</caption><thead><tr><th rowspan="2">[η<sub>h</sub><sup>inf</sup>, η<sub>h</sub><sup>sup</sup>]</th><th rowspan="2">Theoretical mean</th><th rowspan="2">Mean of the 10 samples of η<sub>h</sub></th><th colspan="2">θ<sub>2</sub> = 1 %</th><th colspan="2">θ<sub>2</sub> = 0.001 %</th></tr><tr><th>E[PHSSR] (%)</th><th>Error on E[PHSSR] for a confidence level = 99% (%)</th><th>E[PHSSR] (%)</th><th>Error on E[PHSSR] for a confidence level = 99% (%)</th></tr></thead><tbody><tr><td>[0.0117,0.0143]</td><td>0.0130</td><td>0.0129</td><td>48.82</td><td>0.2833</td><td>48.37</td><td>0.1423</td></tr><tr><td>[0.0143,0.0175]</td><td>0.0159</td><td>0.0161</td><td>49.84</td><td>0.1776</td><td>49.32</td><td>0.1091</td></tr><tr><td>[0.0175,0.0214]</td><td>0.0194</td><td>0.0193</td><td>50.07</td><td>0.0994</td><td>49.85</td><td>0.0224</td></tr><tr><td>[0.0214,0.0261]</td><td>0.0237</td><td>0.0234</td><td>50.78</td><td>0</td><td>50.07</td><td>0.0005</td></tr><tr><td>[0.0261,0.0319]</td><td>0.0290</td><td>0.0284</td><td>50.15</td><td>0.1765</td><td>49.90</td><td>0.0280</td></tr></tbody></table>

Fig. F.12 displays the evolutions of E[PHSSR], PHSSR<sup>max</sup> and PHSSR<sup>min</sup> as functions of $\hat{\eta}_h$. In this case, the maximum value of the objective function is 50.07%, which corresponds to an optimal theoretical value of $\hat{\eta}_h$ equal to 0.0237. The number of simulations is equal to $5 \times 10 \times 10000 = 500000$.

<figure><img src="image_126.png" /><figcaption>**FIGURE F.12.** Evolutions of E[PHSSR], PHSSR<sup>max</sup> and PHSSR<sup>min</sup> as functions of $\hat{\eta}_h$ obtained with the reference method, the ME-gPC based method with $\theta_2 = 1\%$, and the ME-gPC based method with $\theta_2 = 0.001\%$. The black segments denote the limits of the intervals. The reference results are obtained as described in Sect. F.7.1.1.</figcaption></figure>

### F.7.1.2 Optimization with the ME-gPC method

Tab. F.4 shows the results of the optimization of $\hat{\eta}_h$ by the ME-gPC based method. The gPC order is again chosen as $p = 1$, $\theta_1$ is set to $1 \cdot 10^{-3}$ and two values of $\theta_2$ are tested: $\theta_2 = 1\%$ and $\theta_2 = 0.001\%$. The design space is the same as for the reference optimization and the same quantities are considered.
Tab. F.5 compares the results obtained with the three strategies (reference optimization and ME-gPC based optimizations with $\theta_2 = 1\%$ and $\theta_2 = 0.001\%$). The corresponding evolutions of $\mathbb{E}[\text{PHSSR}]$, $\text{PHSSR}^{\text{max}}$ and $\text{PHSSR}^{\text{min}}$ according to $\hat{\eta}_h$ are displayed in Fig. F.12.

In the case of $\theta_2 = 1\%$, the number of simulations is equal to 4613, leading to a computational cost reduction of 98.98%. The maximum value of the objective function is here 50.78%, which corresponds to the same optimal value $\hat{\eta}_h = 0.0237$ as for the reference case. The relative error of the maximum PHSSR mean value with respect to the reference is $E_{\mathbb{E}[\text{PHSSR}]^{\text{max}}} = 1.41\%$.

For $\theta_2 = 0.001\%$, the number of simulations is increased to 8006, which still induces a high computational cost reduction of 98.38%. The maximum value of the objective function is equal to that of the reference case, 50.07%, leading to a relative error $E_{\mathbb{E}[\text{PHSSR}]^{\text{max}}} = 0\%$. The corresponding optimal value of $\hat{\eta}_h$ is again equal to 0.0237.

As for the optimization with deterministic design variables, it is again shown that a proper choice of the precision threshold $\theta_2$ enables the reduction of the error levels. However, an error reduction leads to an increase of the required simulation number (from 4613 to 8006), and thus the retained precision level depends on the accepted computational cost. In both cases presented above, the computational costs with the ME-gPC optimization remains much weaker than that of the reference method (500 000). Moreover, it can be noticed that the optimal values of the mean of $\eta_h$, $\hat{\eta}_h = 0.0237$, and of the PHSSR mean, $\mathbb{E}(\text{PHSSR}) = 50.07\%$, are the same as those obtained with the deterministic optimization. This result is logical because the curve of the PHSSR as a function of $\eta_h$ varies very little around the optimal value. Thus, the PHSSR mean remains practically constant.

<table><caption>TABLEAU F.5. Comparison between the reference optimization and the ME-gPC based optimization. The design parameter is $\hat{\eta}_h$. The reference results are obtained as described in Sect. F.7.1.1.</caption><thead><tr><th></th><th>Nb of simulations</th><th>Computational time</th><th>Computational cost reduction (%)</th><th>E[PHSSR]<sup>max</sup></th><th>E<sub>E</sub>[PHSSR]<sup>max</sup></th><th>η<sub>h</sub><sup>opt</sup></th></tr></thead><tbody><tr><td>Reference</td><td>500 000</td><td>7.86 days</td><td>-</td><td>50.07% ± 0.02%</td><td>-</td><td>0.0237</td></tr><tr><td>ME-gPC<br>(θ<sub>1</sub> = 1 · 10<sup>-3</sup><br>and θ<sub>2</sub> = 1%)</td><td>4 613</td><td>1.92 hours</td><td>98.98</td><td>50.78% ± 2%</td><td>1.41%</td><td>0.0237</td></tr><tr><td>ME-gPC<br>(θ<sub>1</sub> = 1 · 10<sup>-3</sup><br>and<br>θ<sub>2</sub> = 0.001%)</td><td>8 006</td><td>3.29 hours</td><td>98.38</td><td>50.07% ± 0.02%</td><td>0%</td><td>0.0237</td></tr></tbody></table>

## F.7.2 Optimization of the NES nonlinear stiffness

### F.7.2.1 Reference optimization

The last optimization problem proposed in this paper (still represented by (F.18)) considers a single design variable (hyperparameter) $\hat{\varphi}_h$ which is here the mean of the NES nonlinear stiffness $\varphi_h$; similarly to the previous optimization problem, $\varphi_h$ is a stochastic parameter supposed to follow a uniform law. The NES damping coefficient $\eta_h$ is set to 0.024.

The methodology is similar to that of the previous section with a design space divided into $N_I = 5$ intervals as shown in Tab. F.6. The lower and upper bounds of the intervals are again defined as $\varphi_h^{\text{inf}} = \hat{\varphi}_h - 0.1\hat{\varphi}_h$ and $\varphi_h^{\text{sup}} = \hat{\varphi}_h + 0.1\hat{\varphi}_h$, where the mean value $\hat{\varphi}_h$ is chosen so that the upper bound of one interval coincides with the lower bound of the next interval.

As previously, $N_d = 10$ random samples of stiffness values $\varphi_h$ are drawn in each interval according to the uniform probability law, and for each sample the mitigation limit and the PHSSR are directly determined from the numerical simulation of (F.2) using 10000 values of the friction coefficient $\mu$. The objective function to be maximized is the mean $\mathbb{E}[\text{PHSSR}]$ of the 10 PHSSR values obtained in each interval.

In section F.6.2 it has been shown that the PHSSR is a decreasing function of $\varphi_h$, which implies that the optimal value is the lowest value of $\varphi_h$ for which the jump criterion (F.16) is met. It is therefore possible that a part of the interval where this optimal value lies (in concrete terms here, some of the 10 samples in the first interval [38 438, 46 980] N·kg<sup>-1</sup>·m<sup>-3</sup>) corresponds to $\varphi_h$ values for which the jump criterion (F.16)
is not respected. In that case, 10 new samples of $\phi_h$ values are generated in the part of the interval (here [38 438, 46 980]) for which (F.16) is satisfied.

Tab. F.6 shows the reference results of the $\hat{\phi}_h$ optimization. Similarly to the previous section, the mean of the 10 samples $\phi_h$ is calculated for each interval and compared to the theoretical $\hat{\phi}_h$ value. The mean value $E[\text{PHSSR}]$ is also indicated, along with the error on $E[\text{PHSSR}]$ for a confidence level of 99%. With very low error levels (lower than 0.03%), the number of $\phi_h$ samples $N_d = 10$ in each interval is again found sufficient. This reference optimization has required a simulation number equal to $5 \times 10 \times 10 000 = 500 000$.

<table><caption>**TABLEAU F.6.** Results of the reference optimization (see Sect. F.7.2.1). The design parameter is $\hat{\phi}_h$.</caption><thead><tr><th></th><th>Theoretical<br>mean φ̂<sub>h</sub></th><th>Mean of<br>10 samples<br>of φ<sub>h</sub></th><th>E[PHSSR]<br>(%)</th><th>Error on<br>E[PHSSR] for a<br>confidence level<br>= 99% (%)</th></tr></thead><tbody><tr><td>[38 438, 46 980]</td><td>42 709</td><td>46 678</td><td>52.32</td><td>0.0004</td></tr><tr><td>[46 980, 57 420]</td><td>52 200</td><td>52 992</td><td>51.95</td><td>0.0299</td></tr><tr><td>[57 420, 70 180]</td><td>63 800</td><td>62 746</td><td>51.53</td><td>0.0296</td></tr><tr><td>[70 180, 85 776]</td><td>77 978</td><td>78 942</td><td>51.07</td><td>0.0095</td></tr><tr><td>[85 776, 104 837]</td><td>95 306</td><td>94 876</td><td>50.76</td><td>0.0200</td></tr></tbody></table>

The evolutions of $E[\text{PHSSR}]$, $\text{PHSSR}^{\text{max}}$ and $\text{PHSSR}^{\text{min}}$ according to $\hat{\phi}_h$ are displayed in Fig. F.13. In this case, the maximum value of the objective function is 52.32% which corresponds to an optimal value of $\hat{\phi}_h$ equal to 46678. The very low error on the $E[\text{PHSSR}]$ for a 99% confidence level obtained in this case is due to the fact that the PHSSR values are close in the interval.

<figure><img src="image_127.png" /><figcaption>**FIGURE F.13.** Evolutions of $E[\text{PHSSR}]$, $\text{PHSSR}^{\text{max}}$ and $\text{PHSSR}^{\text{min}}$ as functions of $\hat{\phi}_h$ with the reference method, the ME-gPC based method with $\theta_2 = 1\%$ and the ME-gPC based method with $\theta_2 = 0.001\%$. The black segments denote the limits of the intervals. The reference results are obtained as described in Sect. F.7.2.1.</figcaption></figure>

### F.7.2.2 Optimization with the ME-gPC method

Tab. F.7 shows the results of the optimization of $\hat{\phi}_h$ based on the ME-gPC, with again a gPC order $p = 1$, $\theta_1 = 4.8 \cdot 10^{-3}$ and two values of $\theta_2$: $\theta_2 = 1\%$ and $\theta_2 = 0.001\%$. The design space is discretized as for the reference optimization. For several intervals the error on $E[\text{PHSSR}]$ for a confidence level of 99% is null, meaning that the 10 PHSSR values corresponding to the 10 samples are identical. The evolutions of $E[\text{PHSSR}]$, $\text{PHSSR}^{\text{max}}$ and $\text{PHSSR}^{\text{min}}$ according to $\hat{\phi}_h$ for the cases $\theta_2 = 1\%$ and $\theta_2 = 0.001\%$ are visible in Fig. F.13 along with the reference results.

The comparison of the different strategies to optimize the NES nonlinear stiffness (that is, the reference optimization and the two ME-gPC based optimizations with $\theta_2 = 1\%$ and $\theta_2 = 0.001\%$) is presented in Tab. F.8. In the case of $\theta_1 = 4.8 \cdot 10^{-3}$ and $\theta_2 = 1\%$, the number of simulations is equal to 3248, leading to a computational cost reduction of 99.21%. The maximum value of the objective function is 52.73%, which
<table><caption>TABLEAU F.7. Results of the optimization with the ME-gPC method. The design parameter is $\hat{\varphi}_h$.</caption><thead><tr><th rowspan="2">[φ<sub>h</sub><sup>inf</sup>, φ<sub>h</sub><sup>sup</sup>]</th><th rowspan="2">Theoretical mean φ̂<sub>h</sub></th><th rowspan="2">Mean of the 10 samples of φ<sub>h</sub></th><th colspan="2">θ<sub>2</sub> = 1 %</th><th colspan="2">θ<sub>2</sub> = 0.001 %</th></tr><tr><th>E[PHSSR] (%)</th><th>Error on E[PHSSR] for a confidence level = 99% (%)</th><th>E[PHSSR] (%)</th><th>Error on E[PHSSR] for a confidence level = 99% (%)</th></tr></thead><tbody><tr><td>[38 438, 46 980]</td><td>42 709</td><td>46 678</td><td>52.73</td><td>0.4721</td><td>52.33</td><td>0.0007</td></tr><tr><td>[46 980, 57 420]</td><td>52 200</td><td>52 992</td><td>52.34</td><td>0</td><td>51.91</td><td>0.0148</td></tr><tr><td>[57 420, 70 180]</td><td>63 800</td><td>62 746</td><td>52.10</td><td>0.2230</td><td>51.48</td><td>0.0399</td></tr><tr><td>[70 180, 85 776]</td><td>77 978</td><td>78 942</td><td>51.56</td><td>0</td><td>51.09</td><td>0.0105</td></tr><tr><td>[85 776, 104 837]</td><td>95 306</td><td>94 876</td><td>51.09</td><td>0.2599</td><td>50.76</td><td>0.0206</td></tr></tbody></table>

<table><caption>TABLEAU F.8. Comparison between the reference optimization and the optimization with the ME-gPC method. The design parameter is $\hat{\varphi}_h$. The reference results are obtained as described in Sect. F.7.2.1.</caption><thead><tr><th></th><th>Nb of simulations</th><th>Computational time</th><th>Computational cost reduction (%)</th><th>E[PHSSR]<sup>max</sup></th><th>E<sub>E</sub>[PHSSR]<sup>max</sup></th><th>φ̂<sub>h</sub><sup>opt</sup></th></tr></thead><tbody><tr><td>Reference</td><td>500 000</td><td>7.53 days</td><td>-</td><td>52.32% ± 0.02%</td><td>-</td><td>46 678</td></tr><tr><td>ME-gPC<br>(θ<sub>1</sub> = 4.8 · 10<sup>-3</sup> and<br>θ<sub>2</sub> = 1%)</td><td>3248</td><td>1.42 hours</td><td>99.21</td><td>52.73% ± 1.9%</td><td>1.50%</td><td>46 678</td></tr><tr><td>ME-gPC<br>(θ<sub>1</sub> = 2 · 10<sup>-3</sup> and<br>θ<sub>2</sub> = 0.001%)</td><td>6947</td><td>3.06 hours</td><td>98.30</td><td>52.33% ± 0.019%</td><td>0.09%</td><td>46 678</td></tr></tbody></table>

corresponds to an optimal value of $\hat{\varphi}_h$ equal to 46 678 N·kg<sup>-1</sup>·m<sup>-3</sup>, and the relative error on the maximum value of the PHSSR mean is $E_{E[PHSSR]}^{max} = 1.5\%$.

When $θ_1 = 2 · 10^{-3}$ and $θ_2 = 0.001\%$, the number of simulations is increased up to 6947, inducing a slightly reduced computational cost reduction of 98.30%. The maximum value of the objective function is now 52.33%, the optimal value of $\hat{\varphi}_h$ 46 678 N·kg<sup>-1</sup>·m<sup>-3</sup>, and the relative error on the maximum value of the PHSSR mean $E_{E[PHSSR]}^{max} = 0.09\%$.

The conclusions given for the optimization of the damping coefficient with the ME-gPC based method hold here : the choice of the precision threshold $θ_2$ enables the control of the error levels but impacts the computational costs, which remain however much lower than that with the reference method (500 000 simulations).

## F.8 CONCLUSION

In this paper, the robust optimization of NES used to mitigate limit cycles of a friction system with an uncertain friction coefficient has been studied. The optimization strategy is based on the maximization of the propensity of the system to be in an attenuated regime (PHSSR), which constitutes the objective function of the problem. Two original approaches are proposed, depending on whether the parameters of the NES are considered as deterministic or as uncertain, but with a known probability distribution. In each case, two optimization strategies are developed, that consist in evaluating the objective function by either a deterministic reference optimization or by a method based on the Multi-Element generalized Polynomial Chaos (the so-called ME-gPC method). The comparison between the two methods shows the efficiency of the ME-gPC method to significantly reduce the computational cost while keeping a good accuracy compared to
the reference optimization.

These optimization approaches have been applied to determine the optimal values of the NES parameters to maximize their ability to passively attenuate limit cycle of oscillations. Concerning the NES nonlinear stiffness, the optimal value is the lowest value for which the discontinuity jump is considered as high enough. Concerning the NES damping coefficient, the curve representing the PHSSR as a function of the damping coefficient is nonlinear and presents a maximum corresponding to the optimal NES damping value. However, the variations of the PHSSR near the optimal value are very small, which makes it more difficult (but also less crucial) to precisely locate this optimal value.

The general nature of the proposed strategy suggests that it could be used for the robust optimization of any vibration absorber, providing that the transition from the attenuated state to the non-attenuated state passes through a discontinuity in the amplitude profile of the system as a function of the bifurcation parameter considered.

In future works, several issues should be tackled. First, as one of the difficulties of the proposed approach is to choose relevant values for the thresholds, it could be interesting to develop automatically adaptive thresholds. Another issue lies in the extension of the proposed method to several uncertain parameters. It can also be envisaged to study the efficiency of the method for a model with a higher number of degrees of freedom.
# ANALYTICAL PREDICTION OF DELAYED HOPF BIFURCATIONS IN A SIMPLIFIED STOCHASTIC MODEL OF REED MUSICAL INSTRUMENTS

Baptiste BERGEOT and Christophe VERGEZ

## ABSTRACT

This paper investigates the dynamic behavior of a simplified single reed instrument model subject to a stochastic forcing of white noise type when one of its bifurcation parameters (the dimensionless blowing pressure) increases linearly over time and crosses the Hopf bifurcation point of its trivial equilibrium position. The stochastic slow dynamics of the model is first obtained by means of the stochastic averaging method. The resulting averaged system reduces to a non-autonomous one-dimensional Itô stochastic differential equation governing the time evolution of the mouthpiece pressure amplitude. Under relevant approximations the latter is solved analytically treating separately cases where noise can be ignored and cases where it cannot. From that, two analytical expressions of the bifurcation parameter value for which the mouthpiece pressure amplitude gets its initial value back are deduced. These special values of the bifurcation parameter characterize the effective appearance of sound in the instrument and are called deterministic dynamic bifurcation point if the noise can be neglected and stochastic dynamic bifurcation point otherwise. Finally, for illustration and validation purposes, the analytical results are compared with direct numerical integration of the model in both deterministic and stochastic situations. In each considered case, a good agreement is observed between theoretical results and numerical simulations, which validates the proposed analysis.

**Keywords:** Single reed instruments, Self-sustained oscillations, Dynamic bifurcation, Bifurcation delay, Stochastic averaging

## LIST OF MAIN SYMBOLS

The next list is not exhaustive, it describes only the main symbols that will be later used within the body of the document.

$p_n, p, u$ Physical unknowns in the reed instrument model (n<sup>th</sup> modal component, pressure and volume flow at the input of the instrument respectively)

$x_t, x_y$ Amplitude of p in the amplitude/phase representation (different subscript is used to stress the dependence)

$A(y)$ Intermediate quantity used throughout the paper and defined by Eq. (G.34)

$\gamma, \gamma_t$ Dimensionless blowing pressure (the subscript (t) indicates that $\gamma$ varies with time)

$\hat{\gamma}^{st}$ Value of $\gamma$ corresponding to a static Hopf bifurcation

$y_t$ Time-varying dimensionless blowing overpressure with respect to the static Hopf bifurcation (= $\gamma_t - \hat{\gamma}^{st}$)

$\hat{y}_{det}^{dyn}$ Analytical prediction of the value of y at which the system undergoes a dynamic Hopf bifurcation in the deterministic case

$y_0$ Value of $y_t$ at $t = 0$

$\hat{y}_{stoch,a}^{dyn}$ Analytical prediction of the value of y at which the system undergoes a dynamic Hopf bifurcation in the stochastic case (the subscript (a) means first level of approximation)

$\hat{y}_{stoch,b}^{dyn}$ Analytical prediction of the value of y at which the system undergoes a dynamic Hopf bifurcation in the stochastic case (the subscript (b) means second level of approximation, less accurate than $\hat{y}_{stoch,a}^{dyn}$)

$\nu, \sigma$ Magnitude of the stochastic forcing ($\sigma = \nu/\sqrt{2}$)

$\epsilon$ Rate of linear increasing of the bifurcation parameter

$\alpha_i, \omega_i, F_i$ Modal parameters of the i<sup>th</sup> mode (damping, eigenfrequency, modal factor respectively)
ζ

Embouchure parameter (physical model of clarinet)

G.1 INTRODUCTION

Musical instruments are non-linear dynamical systems. An important specificity is that sound production in a musical context corresponds to a time-varying control. Indeed, control parameters are modified continuously by the instrument player. To give just one example, wind instruments players control air pressure in their mouth with variations over time finely tuned to produce the desired sound effect. On the other hand, when studying the corresponding mathematical models of sound production by musical instruments, the control parameters are considered constant over time most of the time. In this respect, the instruments are modeled by autonomous nonlinear systems of differential equations (ODEs) having, among other solutions, at least one trivial equilibrium position (corresponding to silence) and periodic solutions (corresponding to musical notes). The mouth pressure is then a bifurcation parameter of this ODEs system : when it is small the trivial solution is stable and, beyond a precise value called static bifurcation point (“static” because the bifurcation parameter is constant over time), the trivial solution loses its stability through a Hopf bifurcation giving rise to a periodic solution. In musical acoustics literature, the static bifurcation point of the trivial solution is often called oscillation threshold [49]. In this “static” context, periodic solutions and their stability can be also determined using for example the harmonic balance method [42, 71] or orthogonal collocation [81, 196]. Control parameters are changed over time only in the very special context of time simulations, e.g. when the objective is sound synthesis, not model behavior analysis.

The purpose of this paper is not to question these methods (with constant control parameters) which have shown their relevance (experimentally and numerically) when considering steady-state oscillation regimes (i.e. excluding all forms of transients). On the contrary, it is more a question of shedding light on the appearance of a sound when a time-varying control parameter (more specifically, the mouth pressure) crosses the static (Hopf) bifurcation point, as it is the case in playing situations. Moreover, experiments with a blowing machine and a clarinet-like instrument have shown in this case the existence of a delay at the bifurcation [a4] : when the pressure in the musician’s mouth is increased linearly, the start of oscillations is observed for a pressure value greater than the static bifurcation point. This particular value of the mouth pressure for which the oscillations actually appear is called dynamic bifurcation point.

An analytical study to explain this observation has already been carried out on a discrete-time model of the clarinet<fnref n="1" /> [a5, a6]. Considering a blowing pressure linearly increasing over time, the following results have been obtained in these works, which are typical of what is known in the field of dynamical bifurcation of discrete-time systems. In fact, the behavior of the bifurcation delay depends on the possibility of ignoring or not the presence of noise in the model. We are referring here only to the presence of an additive white noise. Sometimes even a noise with a very low amplitude coming from the rounding errors of the computer in numerical simulations must be taken into account. If the noise can be ignored the bifurcation delay depends only on the initial value of the linearly increasing mouth pressure, the farthest it starts below the static bifurcation point, the larger the delay. If the noise is no longer negligible, then the bifurcation delay loses the dependence on the initial condition and logically becomes dependent on the noise level but also on the increase rate of the mouth pressure. In this case, the bifurcation delay decreases with the noise level and increases with the increase rate of the mouth pressure. These properties of the bifurcation delay shown in these works are also found in other works in applied mathematics and physics whether for discrete or continuous-time systems with an additive white noise [11, 12, 23, 106]. The strong dependence of the bifurcation delay on noise makes essential to take it into account in the modeling a real life system like a musical instrument.

The discrete-time clarinet model considered in our previous cited works is known for its simplicity and its ability to explain some phenomena observed experimentally. However, the simple formulation of this model is achieved under very strict simplifying hypotheses. The most important ones are frequency independent losses, a cylindrical geometry, idealized radiation ...

Another approach, based on continuous-time models including ODE’s, allows an easier consideration of refinements in the equations, leading to models with different levels of complexity [182]. The simplest model consists in retaining only one acoustic mode of the air-column in the instrument while assuming that the cane reed is driven instantaneously by the pressure difference between the mouth of the player and the input of the instrument [37]. The work presented in this paper is based on this model which is here considered with a blowing pressure linearly increasing over time and a stochastic forcing of white noise type.

The questions addressed in this paper are : How does the continuous-time clarinet model behave when the blowing pressure increases linearly, compared to the classical case of a blowing pressure constant over time ?

<footnote n="1">1. In this case the instrument is modeled by a difference equation similar to the logistic map. This model has been thoroughly studied in musical acoustics for control parameters constant over time [48, 194, 195].</footnote>
<figure><img src="image_128.png" /><figcaption>FIGURE G.1. Main steps of the proposed approach as a block diagram.</figcaption></figure>

How do the results change when random fluctuations in the source term of the model are added (which may represent an idealized version of the turbulent noise due to the airflow)? More precisely, through relevant approximations, the objective is to solve the studied model analytically by considering the general framework of the stochastic differential equations in order to predict the dynamic bifurcation points of the model.

The paper is organized as follows. The single reed instrument model is presented in Section G.2. Section G.2.1 gives the equations of motion of the classical deterministic one-mode model of single reed instruments and recalls the expression of its static bifurcation point. Section G.2.2 presents the model with, in addition, a stochastic forcing of white noise type, a linearly increasing blowing pressure and pertinent rescaling. In Section G.3, the stochastic averaging method is used to derive the slow dynamics of the model. Through simplifying assumptions, the slow dynamics is solved in Section G.4. The method is based on treating separately cases where noise can be ignored and cases where it cannot. An analytical expression of the dynamic bifurcation point is obtained in each case. In Section G.5, for illustration and validation purposes, these analytical results are compared with direct numerical integration of the model. Finally, concluding remarks and some perspectives are given in Section G.6.

The main steps of the proposed approach are summarized in Fig. G.1.

## G.2 SINGLE REED INSTRUMENT MODEL WITH A WHITE NOISE FORCING

### G.2.1 Deterministic model of single reed instrument and its static bifurcation point

In this section the single reed instrument model is recalled as it is generally discussed in the literature, i.e. it is deterministic and the control parameters are constant over time.

Sound production by single reed instruments is classically modeled through the nonlinear coupling of two linear sub-systems [14, 37, 66] : the cane reed and the air-column inside the instrument. While blowing
air through the reed channel into the instrument, the musician provides a quasi-static source of energy. The instrument and the player constitute an autonomous dynamical system. When the trivial equilibrium solution of this system becomes unstable, a sound is produced [67, 181, 215].

Since the lowest resonance frequency of the reed is one order of magnitude higher than the sound frequency for many notes, the reed is often modeled as a lossless stiffness spring [10, 155]. Therefore, the position of the reed relative to rest (which determines the opening of the reed channel) is proportional to the pressure drop across the reed, i.e. the pressure difference between the mouth and the mouthpiece of the instrument. The linear pressure response of the air column P to the volume flow U through the reed channel is given in the frequency domain through the input impedance of the air column Z :

$$P(\omega) = Z(\omega)U(\omega), \qquad (G.1)$$

where $\omega$ is the angular frequency. The contribution at the input of the instrument of the (infinite) series of modes of the air column is taken into account in $Z(\omega)$. For computational reasons, the series is truncated to $N$ modes, where $N$ is an integer :

$$Z(\omega) = \sum_{n=1}^{N} F_n \frac{j\omega}{\omega_n^2 + j\omega\alpha_n\omega_n - \omega^2}, \qquad (G.2)$$

with $F_n$, $\omega_n$ and $\alpha_n$ the modal parameters, respectively the modal factor, the resonance angular frequency and the inverse of the quality factor of the $n^{th}$ peak of the impedance (corresponding to the $n^{th}$ mode of the air column). Eq. (H.22) can be written in the time domain :

$$\ddot{p}_n + \alpha_n \omega_n \dot{p}_n + \omega_n^2 p_n = F_n \dot{u}, \quad \forall n \in [1, N], \qquad (G.3)$$

with $u$ the inverse Fourier transform of $U$ and $p_n$ is such that $p = \sum_{i=1}^{N} p_n$, where $p$ is the inverse Fourier transform of $P$ [181] and corresponds to the time evolution of the mouthpiece pressure.

The volume flow through the reed channel is related nonlinearly to the reed channel opening and the pressure difference between the mouth and the mouthpiece [50, 98]. A polynomial expansion of this relation is often written in the neighborhood of the equilibrium solution (i.e. the mean flow) [111] :

$$u(t) = u_{eq} + c_1 p(t) + c_2 p(t)^2 + c_3 p(t)^3, \qquad (G.4)$$

with $u_{eq} = \zeta(1 - \gamma)\sqrt{\gamma}$ the mean volume flow, $c_1 = \zeta \frac{3\gamma - 1}{2\gamma^2}$, $c_2 = -\zeta \frac{3\gamma + 1}{8\gamma^2}$ and $c_3 = -\zeta \frac{\gamma + 1}{16\gamma^2}$, where $\gamma$ is the dimensionless pressure in the mouth of the musician and $\zeta$ a dimensionless parameter accounting for many embouchure parameters. By Eq. (G.4), Eq. (H.23) can be written using only the pressure $p$ as follows

$$\ddot{p}_n + \alpha_n \omega_n \dot{p}_n + \omega_n^2 p_n + F_n \dot{p} f(p, \gamma) = 0, \quad \forall n \in [1, N], \qquad (G.5)$$

where $f(p, \gamma) = -\frac{\partial u}{\partial p}$.

A minimal model of a reed instrument including a single mode of the air-column is obtained by stating $N = 1$. In this case (H.25) becomes

$$\ddot{p} + \alpha_1 \omega_1 \dot{p} + \omega_1^2 p + F_1 \dot{p} f(p, \gamma) = 0 \qquad (G.6)$$

In this case ($N = 1$), since $p_1 = p$, $p_1$ is replaced by $p$ in Eq. (H.26). Note that $p$ and $u$ are dimensionless and $F_1$ unit is $s^{-1}$. This is clearly a minimal yet useful model of sound production in reed instruments. Indeed, it takes into account the two main control parameters adjusted by the musician and describes the physical mechanism through which sound emerges from equilibrium (i.e. silence) when a resonance of the air column is excited by an incoming flow.

In this paper the bifurcation parameter under consideration is $\gamma$. The stability of the trivial equilibrium solution $p = 0$ with respect to $\gamma$ is classically analyzed by looking at the sign of the eigenvalues real parts of the Jacobian matrix of Eq. (H.26) written in the state-space form. This leads to the following expression for the static Hopf bifurcation point

$$\hat{\gamma}^{\text{st}} = \frac{1}{3} + \frac{2\alpha_1\omega_1(\alpha_1\omega_1 + \sqrt{\alpha_1^2\omega_1^2 + 3\zeta^2F_1^2})}{9\zeta^2F_1^2} \qquad (G.7)$$

corresponding to the value of $\gamma$ for which the two complex conjugate eigenvalues become with positive real parts. Note that in the lossless case (i.e. $\alpha_1 = 0$) the static bifurcation parameter is $\hat{\gamma}^{\text{st}} = \frac{1}{3}$. Static bifurcation point is known in the literature of acoustics of musical instruments as the "oscillation threshold" of the instrument. Expressions equivalent to Eq. (G.7) are given for example by Kergomard et al. [112] and Silva et al. [183].

Compared to the model given by (H.26), this article studies the case of a linearly increasing blowing pressure $\gamma$, with additional stochastic excitation as detailed in Sect. G.2.2.
### G.2.2 One-mode stochastic single-reed instrument model with a linearly increasing blowing pressure

Hereafter, the formalism of stochastic differential equations is used in the framework of the Itô stochastic calculus. A good description of these concepts can be found for example in [153].

The one-mode model described by Eq. (H.26) is now considered with a linearly increasing blowing pressure $\gamma = \hat{\epsilon}t + \gamma_0$. The variation of the blowing pressure is assumed to be small during a period $T_1 = \frac{2\pi}{\omega_1}$ which is ensured by $0 < \hat{\epsilon} \ll 1$. Because one assumes also that $\omega_1 \gg 1$ (without loss a generality in musical context) the time variation of the mouth pressure $\gamma$ is neglected in the time derivative of the mean flow (G.4). Moreover, the model is now forced by a stochastic excitation $\hat{\nu}\xi_t$ where $\hat{\nu}$ is the noise level and $\xi_t$ is assumed to be a unitary idealized white noise process with zero mean, i.e.

$$ \mathbb{E}[\xi_t] = 0 \quad \text{and} \quad \mathbb{E}[\xi_t \xi_{t+\tau}] = \delta(\tau) \qquad (G.8) $$

where $\delta$ is the Dirac delta function and $\mathbb{E}[\{\cdot\}]$ denotes the ensemble average<fnref n="2" /> of $\{\cdot\}$. Finally, the one-mode stochastic model is written as follows

$$ \ddot{p}_t + \alpha_1 \omega_1 \dot{p}_t + \omega_1^2 p_t + F_1 \dot{p}_t f(p_t, \gamma_t) = \hat{\nu} \xi_t \qquad (G.9) $$

where the subscript $t$ is used to show the stochastic nature of the differential equation. This is the classical notation in the framework of stochastic differential equations.

From the physical point of view, Eq. (G.9) can be interpreted as a simple model of sound production in reed instruments with turbulent noise contribution taken into account as an additive stochastic noise in the source term.

To facilitate the mathematical developments of the following sections, transformations are performed within Eq. (G.9). First, a new bifurcation parameter $y_t = \gamma_t - \hat{\gamma}^{\text{st}}$ is considered to obtain a system whose trivial solution becomes unstable at $y = 0$ (i.e. $\hat{y}^{\text{st}} = 0$). Then the time rescaling $t \to t' = \omega_1 t$ is introduced. Recalling that a normalized white noise $\xi_t$ is defined as the time derivative of the normalized Wiener process $W_t$ and using the scaling property of Wiener process (see [153], Chap. 2) which states that $W_{\frac{t}{\omega_1}}$ and $\frac{1}{\sqrt{\omega_1}}W_t$ are equivalent, we have

$$ \xi_t = \frac{dW_t}{dt} \sim \omega_1 \frac{dW_{\frac{t'}{\omega_1}}}{dt'} \sim \sqrt{\omega_1} \frac{dW_{t'}}{dt'} = \sqrt{\omega_1} \xi_{t'}. \qquad (G.10) $$

Therefore, by using also $\{\cdot\}$ for the derivation with respect to $t'$ and denoting $t'$ by $t$ for the sake of conciseness, Eq. (G.9) takes the form of the following stochastically forced self-excited oscillator

$$ \ddot{p}_t + h(p_t, \dot{p}_t, y_t) + p_t = \nu \xi_t \qquad (G.11) $$

where

$$ h(p_t, \dot{p}_t, y_t) = \alpha_1 \dot{p}_t + \frac{F_1}{\omega_1} \dot{p}_t f(p_t, y_t + \hat{\gamma}^{\text{st}}) \qquad (G.12) $$

and $\nu = \frac{\dot{\hat{\nu}}}{\omega_1^{3/2}}$. The time evolution of $y_t$ is given by

$$ y_t = \epsilon t + y_0 \qquad (G.13) $$

where $y_0 = \gamma_0 - \hat{\gamma}^{\text{st}}$ and $\epsilon = \frac{\hat{\epsilon}}{\omega_1}$ with $0 < \hat{\epsilon} \ll 1$ (because $\omega_1 \gg 1$). That means that if $y_t$ is a slow variable for the time scale $t$ it is also slow for the faster time scale $\omega_1 t$.

## G.3 EQUATIONS GOVERNING THE STOCHASTIC SLOW DYNAMICS

Following Roberts and Spanos [169] the stochastic averaging method [114, 191], whose general formulation is recalled in Appendix G.8, is used to obtain the stochastic slow dynamics of Eq. (H.27). To apply the method, the response process $(p_t, \dot{p}_t)$ needs to be transformed into a pair of slowly varying processes. To achieve that, as in classical averaging methods of deterministic systems, an amplitude-phase representation is used imposing

$$ p_t = x_t \cos(t + \varphi_t) \qquad (G.14a) $$

$$ \dot{p}_t = -x_t \sin(t + \varphi_t) \qquad (G.14b) $$

<footnote n="2">2. We recall that the ensemble average consists in repeating the same measurement many times, and in calculating the average over them.</footnote>
The desired form of Eq. (G.14b) requires that

$$ \dot{x}_t \cos \phi_t - x_t \dot{\phi}_t \sin \phi_t = 0, \qquad (G.15) $$

with $\phi_t = t + \varphi_t$ which yields

$$ \dot{\varphi}_t = \frac{\dot{x}_t \cos \phi_t}{x_t \sin \phi_t} \quad \text{and} \quad \dot{x}_t = x_t \dot{\varphi}_t \frac{\sin \phi_t}{\cos \phi_t}. \qquad (G.16) $$

Then differentiation of Eq. (G.14b) leads to

$$ \ddot{p}_t = -x_t \cos \phi_t - \dot{x}_t \sin \phi_t - x_t \dot{\varphi} \cos \phi_t. \qquad (G.17) $$

Finally, the substitution of (G.14) and (G.17) into (H.27) and the use of (G.16) yields

$$ \dot{x}_t = h(x_t \cos \phi_t, -x_t \sin \phi_t, y_t) \sin \phi_t - \nu \xi_t \sin \phi_t \qquad (G.18a) $$

$$ \dot{\varphi}_t = h(x_t \cos \phi_t, -x_t \sin \phi_t, y_t) \frac{\cos \phi_t}{x_t} - \nu \xi_t \frac{\cos \phi_t}{x_t} \qquad (G.18b) $$

which has a similar form as (G.70) in Appendix G.8 with

$$ \begin{aligned} \mathbf{x}_t &= (x_t, \varphi_t)^T, \qquad \eta_t = (\xi_t, \xi_t)^T, \\ \mathbf{f}(\mathbf{x}_t, t) &= \begin{pmatrix} h(x_t \cos \phi_t, -x_t \sin \phi_t, y_t) \sin \phi_t \\ \frac{1}{x_t} h(x_t \cos \phi_t, -x_t \sin \phi_t, y_t) \cos \phi_t \end{pmatrix}, \qquad (G.19) \\ \mathbf{g}(\mathbf{x}_t, t) &= \begin{pmatrix} -\nu \sin \phi_t & 0 \\ 0 & -\nu \frac{\cos \phi_t}{x_t} \end{pmatrix} \end{aligned} $$

with $()^T$ the transpose operator.

Following the approach described in Appendix G.8, the drift vector **m** (see Eq. (G.72)) is first computed. The first term $T^{\text{av}}\{\mathbf{f}\} = (F(x_t, y_t), G(x_t, y_t))^T$, with

$$ F(x_t, y_t) = \left( \frac{F_1}{\omega_1} \zeta \frac{3(y_t + \hat{\gamma}^{\text{st}}) - 1}{4(y_t + \hat{\gamma}^{\text{st}})^{1/2}} - \frac{\alpha_1}{2} \right) x_t + \frac{F_1}{\omega_1} \zeta \frac{3(\hat{\gamma}^{\text{st}} + y_t + 1)}{128 (\hat{\gamma}^{\text{st}} + y_t)^{5/2}} x_t^3 \qquad (G.20) $$

$$ G(x_t, y_t) = 0, \qquad (G.21) $$

corresponds to the classical (deterministic) Bogoliubov-Krylov averaging of the vector function **f**, i.e. the time average is performed assuming that $x_t$ and $\phi_t$ are slow variables with respect to the unitary eigenfrequency of the considered dimensionless simplified reed instrument model (H.27) which therefore corresponds to a period equal to $2\pi$.

Then, the second part of the drift vector **m** (see again Eq. (G.72)) is determined as

$$ \begin{aligned} T^{\text{av}} & \left\{ \int_{-\infty}^{0} \mathbb{E} \left[ \left( \frac{\partial(\mathbf{g}\eta)}{\partial \mathbf{x}} \right)_t (\mathbf{g}\eta)_{t+\tau} \right] d\tau \right\} \\ &= \frac{1}{2\pi} \int_{0}^{2\pi} \int_{-\infty}^{0} \mathbb{E} \left[ \left( \frac{\partial(\mathbf{g}\eta)}{\partial \mathbf{x}} \right)_t (\mathbf{g}\eta)_{t+\tau} \right] d\tau dt \\ &= \frac{1}{2\pi} \int_{-\infty}^{0} \int_{0}^{2\pi} \delta(\tau) \left( \frac{\nu^2 \cos \phi_t \cos \phi_{t+\tau}}{-\nu^2 \sin \phi_t \cos \phi_{t+\tau} + \cos \phi_t \sin \phi_{t+\tau}} \frac{x_t}{x_t^2} \right) dt d\tau \\ &= \left( \frac{\nu^2}{4x_t} \right) \end{aligned} \qquad (G.22) $$

in which the symbol $\mathbb{E}$ disappeared due to definitions (G.8) and one used the fact that $\int_{-\infty}^{0} \delta(\tau) \cos(\tau) d\tau = \frac{1}{2} \int_{-\infty}^{+\infty} \delta(\tau) \cos(\tau) d\tau = \frac{1}{2}$. The final expression of the drift vector is therefore

$$ \mathbf{m} = \begin{pmatrix} F(x_t, y_t) + \frac{\nu^2}{4x_t} \\ 0 \end{pmatrix}. \qquad (G.23) $$
Now the expression of the diffusion matrix $\boldsymbol{\sigma}$ is determined from (G.73) starting by computing

$$ 
\begin{aligned}
T^{\text{av}} & \left\{ \int_{-\infty}^{+\infty} \mathbb{E} \left[ (\mathbf{g}\boldsymbol{\eta})_t (\mathbf{g}\boldsymbol{\eta})_{t+\tau}^T \right] d\tau \right\} \\
&= \frac{1}{2\pi} \int_0^{2\pi} \int_{-\infty}^{+\infty} \mathbb{E} \left[ (\mathbf{g}\boldsymbol{\eta})_t (\mathbf{g}\boldsymbol{\eta})_{t+\tau}^T \right] d\tau dt \\
&= \int_{-\infty}^{+\infty} \left( \frac{\frac{\nu^2 \delta(\tau) \cos(\tau)}{2}}{\frac{\nu^2 \delta(\tau) \sin(\tau)}{2x_t}} - \frac{\frac{\nu^2 \delta(\tau) \sin(\tau)}{2x_t}}{\frac{\nu^2 \delta(\tau) \cos(\tau)}{2x_t^2}} \right) d\tau \\
&= \begin{pmatrix} \frac{\nu^2}{2} & 0 \\ 0 & \frac{\nu^2}{2x_t^2} \end{pmatrix}.
\end{aligned} 
\quad (G.24) $$

in which again definitions (G.8) have been used. From (G.24) a possible solution of (G.73) is

$$ \boldsymbol{\sigma} = \begin{pmatrix} \frac{\nu}{\sqrt{2}} & 0 \\ 0 & \frac{1}{x_t} \frac{\nu}{\sqrt{2}} \end{pmatrix}. \quad (G.25) $$

Therefore, according to (G.71) the equations governing the stochastic slow dynamics of (H.27) are

$$ dx_t = \left( F(x_t, y_t) + \frac{\nu^2}{4x_t} \right) dt + \frac{\nu}{\sqrt{2}} dW_t \quad (G.26a) $$

$$ d\varphi_t = \frac{1}{x_t} \frac{\nu}{\sqrt{2}} dW_t. \quad (G.26b) $$

in which the time evolution of the amplitude $x_t$ is uncoupled from that of the phase $\varphi_t$.

Assuming a weak noise level (i.e. $\nu \ll 1$) and therefore neglecting the term $\frac{\nu^2}{4x_t}$ in (G.26a) and denoting $\sigma = \frac{\nu}{\sqrt{2}}$ we obtain the following Itô stochastic differential equation for the amplitude of the mouth pressure $p_t$

$$ dx_t = F(x_t, y_t)dt + \sigma dW_t. \quad (G.27) $$

The scaling and time shift properties of Wiener process (see [153], Chap. 2) allow us to write that $\frac{1}{\sqrt{\epsilon}}(W_{\epsilon t+y_0} - W_{y_0}) \sim \frac{1}{\sqrt{\epsilon}}W_{\epsilon t} \sim W_t$ and therefore

$$ \frac{1}{\sqrt{\epsilon}} dW_y \sim dW_t. \quad (G.28) $$

Finally, from (G.27) and (G.28) we obtain the final form of the slow dynamics as

$$ dx_y = \frac{1}{\epsilon} F(x_y, y) dy + \frac{\sigma}{\sqrt{\epsilon}} dW_y. \quad (G.29) $$

where $x_t$ and $y_t$ are now denoted $x_y$ and $y$ respectively to emphasize that Eq. (G.29) is a non-autonomous one-dimensional Itô stochastic differential equation with respect to $x_y$ depending on $y$ which acts as the time variable.

## G.4 ANALYTICAL EXPRESSION OF THE DYNAMIC PITCHFORK BIFURCATION POINTS OF THE AVERAGED SYSTEM

The supercritical Hopf bifurcation at $y = 0$ for the initial (non averaged) system (H.27) becomes a supercritical pitchfork bifurcation for Eq. (G.27). Therefore, if the averaged system (G.27) is a good approximation of (H.27), the dynamic Hopf bifurcation point of Eq. (H.27) is very close to the dynamic pitchfork bifurcation point of (G.27). The analytical expression of the latter is obtained in this section in both deterministic and stochastic cases.

As stated by Stocks et al. [190], the dynamic bifurcation point<fnref n="3" /> has two main definitions: it is (1) the value of the bifurcation parameter $y$ when $\sqrt{\mathbb{E}[x_y^2]}$ crosses a predefined threshold $x_{th}^2$ (here the initial value is

<footnote n="3">3. Sometimes called *exit value* in the literature.</footnote>
chosen) or (2) the value of the bifurcation parameter y at which, on average, the random variable $x_y^2$ crosses the same given threshold $x_{th}^2$. In general, these two averaging procedures will lead to two slightly different results. In the paper the first definition is chosen and is clearly recalled below. Other definitions exist (see e.g. [11] and [a6] for a clarinet-like system).

**DEFINITION G.1 (DYNAMIC PITCHFORK BIFURCATION POINT).** In the case of a pitchfork bifurcation and in a deterministic framework, the dynamic bifurcation point is defined as the value of y for which $x_y$ exceeds its initial value $x_{y0}$ (in absolute value). In a stochastic framework, as in this work, the expected value of the squared amplitude is considered and the bifurcation point, denoted $\hat{y}^{dyn}$, is defined as the value of $y_t$ for which $\sqrt{\mathbb{E}[x_y^2]}$ exceeds the initial value $x_{y0}$. The dynamic bifurcation point is larger than the static bifurcation and the difference between them is called **bifurcation delay**.

### G.4.1 Solution of the linearized averaged system

The linearized version of (G.29) with respect to $x_y$ around 0 is considered

$$dx_y = \frac{1}{\epsilon}a(y)x_ydy + \frac{\sigma}{\sqrt{\epsilon}}dW_y \quad (G.30)$$

where, from Eq. (G.20), one has

$$a(y) = \frac{\partial F}{\partial x_y}(0, y) = \frac{F_1}{\omega_1} \zeta \frac{3(y + \hat{\gamma}^{\text{st}}) - 1}{4(y + \hat{\gamma}^{\text{st}})^{1/2}} - \frac{\alpha_1}{2}. \quad (G.31)$$

In order to solve (G.30), the deterministic differential equation associated to the latter is first considered as

$$\frac{dx_y}{dy} = \frac{1}{\epsilon}a(y)x_y \quad (G.32)$$

whose solution, denoted $x^{\det}(y)$, is

$$x^{\det}(y) = x_{y0} e^{\frac{1}{\epsilon}(A(y) - A(y_0))} \quad (G.33)$$

where

$$A(y) = \frac{\zeta F_1(\hat{\gamma}^{\text{st}} + y - 1)\sqrt{\hat{\gamma}^{\text{st}} + y}}{2\omega_1} - \frac{\alpha_1}{2}y. \quad (G.34)$$

is the antiderivative of $a(y)$.

Following a well-known method for solving stochastic differential equations, we apply the Itô formula (see Eq. (G.69) in Appendix G.7) to the function $f(x_y, y) = x_y e^{-\frac{1}{\epsilon}(A(y) - A(y_0))}$ which is constant for the deterministic Eq. (G.32) (indeed through Eq. (G.33) one has $f(x^{\det}(y), y) = x_{y0}$). Using the differential of the product rule, and noting that in the Itô formula the term $\frac{\partial f}{\partial y} + \frac{1}{\epsilon}a(y)x_y \frac{\partial f}{\partial x}$ vanishes, we obtain

$$df(x_y, y) = \frac{\sigma}{\sqrt{\epsilon}}e^{-\frac{1}{\epsilon}(A(y) - A(y_0))}dW_y \quad (G.35)$$

and integrating (G.35) from $y_0$ to y yields the following solution of Eq. (G.30)

$$x_y = x^{\det}(y) + \frac{\sigma}{\sqrt{\epsilon}} e^{\frac{1}{\epsilon}A(y)} \int_{y_0}^{y} e^{-\frac{1}{\epsilon}A(y')} dW_{y'}. \quad (G.36)$$

Note that the first and second terms of the right-hand side of (G.36) are the deterministic and stochastic parts of $x_y$ respectively. The integral in the second term is called an Itô integral and, for a given function g, has the following property $\mathbb{E}[\int g dW_y] = 0$ assuming some properties for the function g (see [153], Chap. 3) that are respected by $e^{-\frac{1}{\epsilon}A(y)}$. Therefore we have

$$\mathbb{E} \left[ \int_{y_0}^{y} e^{-\frac{1}{\epsilon} A(y)} dW_y \right] = 0. \quad (G.37)$$

That means that the expected value of $x_y$ is the solution of the associated deterministic equation (G.32), i.e. $\mathbb{E}[x_y] = x^{\det}(y)$.

Now $\mathbb{E}[x_y^2]$ is computed from (G.36), that yields

$$\mathbb{E}[x_y^2] = \mathbb{E}[(x^{\det}(y))^2]$$
$$ + \mathbb{E} \left[ \frac{\sigma^2}{\epsilon} e^{\frac{2}{\epsilon} A(y)} \left( \int_{y_0}^{y} e^{-\frac{1}{\epsilon} A(y')} dW_{y'} \right)^2 \right] + \mathbb{E} \left[ 2 x^{\text{det}}(y) \frac{\sigma}{\sqrt{\epsilon}} e^{\frac{1}{\epsilon} A(y)} \int_{y_0}^{y} e^{-\frac{1}{\epsilon} A(y')} dW_{y'} \right]. \quad (G.38) $$

Using again the property of the Itô integral which states that $\mathbb{E}[\int gdW_y] = 0$, the third term in the right-hand side of (G.38) is equal to zero. The second term is simplified using the Itô isometry (see Corollary 3.1.7 in [153], Chap. 3) which states that $\mathbb{E}[(\int gdW_y)^2] = \mathbb{E}[\int g^2dy]$ assuming again some properties for the function $g$ that are resected by $e^{-\frac{1}{\epsilon}A(y)}$. Therefore we obtain the final expression of the expected value of the squared amplitude

$$ \mathbb{E}[x_y^2] = \mathcal{D}(y) + \mathcal{S}(y) \quad (G.39) $$

where

$$ \mathcal{D}(y) = (x^{\text{det}}(y))^2 \quad (G.40a) $$

$$ \mathcal{S}(y) = \frac{\sigma^2}{\epsilon} e^{\frac{2}{\epsilon} A(y)} \int_{y_0}^{y} e^{-\frac{2}{\epsilon} A(y)} dy \quad (G.40b) $$

are respectively the deterministic and the stochastic parts of the expected value of the squared amplitude.

From (G.39) and according to [19], three regimes can be distinguished for a system such as (G.29) under-going a dynamic pitchfork bifurcation :

* Regime I : $S(y) \ll \mathcal{D}(y)$. In this case the noise level is so small that it can be neglected and the problem is identical to the deterministic case (G.32) which undergoes the greatest possible bifurcation delay.

* Regime II : $S(y) \gg \mathcal{D}(y)$ with a noise level not too high. The noise can no longer be neglected but it remains small enough for linearized model (G.30) to remain valid. A bifurcation delay still exists but it is reduced compared to Regime I.

* Regime III : $S(y) \gg \mathcal{D}(y)$ with a high noise level. In this situation the behavior of the system is dominated by noise and the trajectory of $x_y$ leaves the neighborhood of zero before the static bifurcation point $y = 0$ is reached. In such a case the linear approximation is therefore not valid anymore and, as highlighted by Berglung and Gentz [19], the notion of bifurcation delay becomes meaningless.

The domains of existence of each regime are given explicitly in Sect. G.4.4.

In the two following parts, we concentrate on Regimes I and II, giving the analytical expression of the dynamic bifurcation points considering each regime separately.

## G.4.2 Dynamic bifurcation point in the deterministic case

In this section we assume that the system is in Regime I and therefore $\mathbb{E}[x_y^2] = \mathcal{D}(y) = (x^{\text{det}}(y))^2$ (with $x^{\text{det}}(y)$ given by Eq. (G.33)). Consequently, from Definition G.1, the deterministic dynamic bifurcation point, denoted $\hat{y}_{\text{det}}^{\text{dyn}}$, is a solution of $x_{y_0} = x_{y_0} e^{\frac{1}{\epsilon}(A(y)-A(y_0))}$ and therefore of

$$ A(y) = A(y_0). \quad (G.41) $$

From (G.34), Eq. (G.41) is solved (details are given in Appendix G.9). It is shown that, in addition to the trivial solution $y = y_0$, Eq. (G.41) has two other solutions. The one that corresponds to the deterministic dynamic bifurcation point is

$$ \hat{y}_{\text{det}}^{\text{dyn}} = X_2^2 - \hat{\gamma}^{\text{st}}, \quad (G.42) $$

where $X_2$ is given by Eq. (G.78). We will see in Fig. G.2 that the maximum of $\hat{y}_{\text{det}}^{\text{dyn}}$ is obtained for $y_0 = -\hat{\gamma}^{\text{st}}$ and in this case Eq. (G.42) reduces to

$$ \hat{y}_{\text{det}}^{\text{dyn}} = \frac{\left( \alpha_1 \omega_1 + \sqrt{\alpha_1^2 \omega_1^2 + 4 \zeta^2 F_1^2} \right)^2}{4 \zeta^2 F_1^2} - \hat{\gamma}^{\text{st}} \quad (G.43) $$

It is important to be aware that even if $\hat{y}_{\text{det}}^{\text{dyn}}$ has a finite value for $y_0 = -\hat{\gamma}^{\text{st}}$, the function $F(x, y)$ in Eq. (G.29) diverge at $y = -\hat{\gamma}^{\text{st}}$. Therefore, Eq. (G.43) must be understood as a limit value for $y_0$ very close to $-\hat{\gamma}^{\text{st}}$ (from above). However, for values of $y$ very close to $-\hat{\gamma}^{\text{st}}$, the linearization of Eq. (G.29) with respect to $x$ (leading to Eq. (G.30)) is valid if $x$ is very close to zero. In other words, Eq. (G.43) is valid to predict the deterministic bifurcation point of a numerical simulation of Eq. (G.29) only if $y_0$ and $x_{y_0}$ are very close to $-\hat{\gamma}^{\text{st}}$ and 0 respectively. This corresponds to the situation where the musician begins to blow gently into the instrument.
<figure><img src="image_129.png" /><figcaption>FIGURE G.2. Deterministic dynamic bifurcation point $\hat{y}_{det}^{dyn}$, given by Eq. (G.42), as a function of the initial condition $y_0$. The opposite of the static bifurcation point $\hat{\gamma}^{st}$, given by Eq. (G.7), is depicted by a vertical dashed line. The red points correspond to the initial conditions used in Fig. G.3. The set of parameters (G.45) is used.</figcaption></figure>

<figure><img src="image_130.png" /><figcaption>FIGURE G.3. Numerical simulations of the deterministic differential equation associated to Eq. (G.29) using a logarithm scale for the vertical axis. The set of parameters (G.45) is used with, in addition, $x(y_0) = 0.01$ (the latter is depicted by a horizontal dashed line) and $y_0 = -0.36, -0.32, \dots, -0.16$. The static bifurcation point $\hat{y}^{st} = 0$ is depicted by a vertical gray dashed line. The dynamic bifurcation point $\hat{y}_{det}^{dyn}$ is defined for a given initial condition by the value of $y$ for which $|x(y)| = |x(y_0)|$, which corresponds graphically to the intersection between the blue curves and the horizontal dashed line. The green and red points are used respectively to highlight the considered initial conditions and the corresponding dynamic bifurcation points predicted by Eq. (G.42).</figcaption></figure>

In the lossless case, i.e. $\alpha_1 = 0$ and $\hat{\gamma}^{st} = \frac{1}{3}$, Eq. (G.42) becomes

$$ \hat{y}_{det}^{dyn} = \frac{1}{2} \left( 1 - y_0 - \sqrt{1 + y_0(2 - 3y_0)} \right) \qquad (G.44) $$

which depends only on the initial value $y_0$. Note first that Eq. (G.44) yields $\hat{y}_{det}^{dyn} = 2/3$ if $y_0 = -\hat{\gamma}^{st} = -\frac{1}{3}$.

This is the largest value of the dynamic bifurcation point. From a physical point of view this means that the oscillations emerge when the blowing pressure $\gamma = \hat{\gamma}^{st} + \hat{y}_{det}^{dyn} = 1$. However this value is known to be the limit of validity of the clarinet model since the reed channel becomes completely closed (an effect not taken into account in the model considered in the work). Therefore, in the lossless case a linear increase of the blowing pressure from an arbitrary small value gives a scenario where the clarinet never plays. This remark can be extended, as shown below, to situations in which damping is taken into account (i.e. $\alpha_1 > 0$).

Then, in cases with and without damping, $\hat{y}_{det}^{dyn}$ does not depend on the slope $\epsilon$ which may seem counte-rintuitive.

The deterministic dynamic bifurcation point $\hat{y}_{det}^{dyn}$, given by Eq. (G.42), is plotted in Fig. G.2 as a function
<figure><img src="image_131.png" /><figcaption>FIGURE G.4. The deterministic bifurcation point $\hat{y}_{det}^{dyn}$ given by Eq. (G.42) as a function of the initial condition $y_0$ for five values of the damping coefficient, i.e. $\alpha_1 = 0, 0.025, \dots, 0.1$ with $\omega_1 = 1000$ rad·s⁻¹, $F_1 = 1200$ s⁻¹ and $\zeta = 0.2$. Red dashed lines indicate the value of $-\hat{\gamma}^{st}$ for each value of $\alpha_n$.</figcaption></figure>

of the initial condition $y_0$ with the following set of parameters :

$$ \epsilon = 0.002, \quad \omega_1 = 1000 \text{ rad·s}^{-1}, \quad \alpha_1 = 0.02, \tag{G.45} \\ F_1 = 1200 \text{ s}^{-1} \quad \text{and} \quad \zeta = 0.2. $$

The figure shows that $\hat{y}_{det}^{dyn}$ is a decreasing function with respect to $y_0$ and, as previously mentioned, the maximum is obtained for $y_0 = -\hat{\gamma}^{st}$. Moreover, the red points correspond to the initial conditions used in Fig. G.3. In the latter, numerical simulations of the deterministic differential deterministic differential equation associated to Eq. (G.29) are shown with a logarithmic scale for the vertical axis. The same parameters are as in Fig. G.2 and the initial values are : $x(y_0) = 0.01$ and $y_0 = -0.36, -0.32, \dots, -0.16$. The dynamic bifurcation point is defined for a given initial condition by the value of y for which $|x(y)| = |x(y_0)|$, which corresponds graphically to the intersection between the blue curves and the horizontal dashed line. The minimum is reached at the static bifurcation point $\hat{y}^{st} = 0$ (depicted by a vertical dashed gray line). One can see that the smaller $y_0$, the smaller the minimum of the trajectory $x(y)$ too and consequently the longer the path to travel before reaching $x(y_0)$ leading to larger values of the deterministic dynamic bifurcation point $\hat{y}_{det}^{dyn}$. The green and red points are used respectively to highlight the considered initial conditions and the corresponding dynamic bifurcation points predicted by Eq. (G.42). This shows that the theoretical results presented in Fig. G.2 predict the deterministic bifurcation points measured on numerical simulations.

In the lossless case the expression of the deterministic dynamic bifurcation point $\hat{y}_{det}^{dyn}$ is given by (G.44) which depends only on the initial value $y_0$. The deterministic bifurcation point $\hat{y}_{det}^{dyn}$ given by Eq. (G.42) is plotted as a function of the initial condition $y_0$ for five values of the damping coefficient in Fig. G.4 and as a function of the parameter $\zeta$ for four values of the damping coefficient $\alpha_1$ and for three values of the initial condition $y_0$ in Fig. G.5. The figures show that in the case with a small damping, the initial value $y_0$ remains the most influential parameter except for the smallest values of the control parameter $\zeta$ and for initial conditions $y_0$ close to the opposite of the static bifurcation point $\hat{\gamma}^{st}$ (see Fig. G.5). In general one has $0.1 < \zeta < 0.4$ for a clarinet, $0.25 < \zeta < 1$ for a saxophone and more for double-reed instruments.

As said previously, the deterministic dynamic bifurcation point corresponds to the largest possible bifurcation delay which holds when the noise can be neglected. In general, because the noise prevents the amplitude $x_y$ to have a very small value when the parameter $y$ reaches its static bifurcation value, the effect of an additive noise is to reduce the delay and to make it lose its dependence on the initial condition $y_0$. The influence of noise is studied in the next section.
<figure><img src="image_132.png" /><figcaption>FIGURE G.5. The deterministic bifurcation point $\hat{y}_{\text{det}}^{\text{dyn}}$ given by Eq. (G.42) as a function of the parameter $\zeta$ for four of the damping coefficient ($\alpha_1 = 0, 0.01, 0.02$ and $0.05$), for three value of the initial condition ($y_0 = -\hat{\gamma}^{\text{st}}, -3\hat{\gamma}^{\text{st}}/4$ and $-\hat{\gamma}^{\text{st}}/2$) with $\omega_1 = 1000 \text{ rad}\cdot\text{s}^{-1}$ and $F_1 = 1200 \text{ s}^{-1}$.</figcaption></figure>

### G.4.3 Dynamic bifurcation point in the stochastic case

Now the system is assumed to evolve according to Regime II, i.e. $E[x_y^2] = S(y)$. In this case, again through Definition G.1 and using Eq. (G.40b), the equation to solve in order to find the dynamic bifurcation point is

$$x_{y_0}^2 = \frac{\sigma^2}{\epsilon} e^{-\frac{2}{\epsilon} A(y)} \int_{y_0}^{y} e^{-\frac{2}{\epsilon} A(y)} dy. \qquad (G.46)$$

The first step is to obtain the approximate expression of the integral

$$I(y_0, y) = \int_{y_0}^{y} e^{-\frac{2}{\epsilon} A(y)} dy. \qquad (G.47)$$

For this purpose $A(y)$ (given by Eq. (G.34)) is expanded in a second order Taylor series around 0 (the static bifurcation point), that leads to

$$A(y) \approx A(0) + y \frac{dA}{dy}(0) + \frac{y^2}{2} \frac{d^2 A}{dy^2}(0) \qquad (G.48a)$$

$$= A(0) + ya(0) + \frac{y^2}{2} a'(0) \qquad (G.48b)$$

where the second term in the right-hand side of (G.48b) vanishes by definition<fnref n="4" /> and $a'(0) = \frac{da}{dy}(0)$. The approximate expression of $I(y_0, y)$ is therefore given at order 2 by

$$\begin{aligned} I(y_0, y) &= \int_{y_0}^{y} e^{-\frac{2}{\epsilon} \left(A(0) + \frac{y^2}{2} a'(0)\right)} dy \\ &= \frac{1}{2} \sqrt{\frac{\pi \epsilon}{a'(0)}} e^{-\frac{2}{\epsilon} A(0)} \left( \operatorname{erf} \left( y \sqrt{\frac{a'(0)}{\epsilon}} \right) - \operatorname{erf} \left( y_0 \sqrt{\frac{a'(0)}{\epsilon}} \right) \right) \end{aligned} \qquad (G.49)$$

where erf is the error function.

<footnote n="4">4. Indeed, at the bifurcation one has $a(0) = \frac{\partial F}{\partial x}(0, 0) = 0$.</footnote>
Since we consider $0 < \epsilon \ll 1$ and because $a'(0)$ has a finite value, one has for an initial value chosen smaller than the static bifurcation point, i.e. for $y_0 < 0$

$$y_0 \sqrt{\frac{a'(0)}{\epsilon}} \ll -1 \Rightarrow \operatorname{erf} \left( y_0 \sqrt{\frac{a'(0)}{\epsilon}} \right) \approx -1 \qquad (G.50)$$

and for $y > 0$ (we are interested in the dynamic bifurcation point which is by definition larger than the static bifurcation point)

$$y \sqrt{\frac{a'(0)}{\epsilon}} \gg 1 \Rightarrow \operatorname{erf} \left( y \sqrt{\frac{a'(0)}{\epsilon}} \right) \approx 1. \qquad (G.51)$$

Therefore, from Eqs. (G.50) and (G.51), Eq. (G.49) reduces to

$$I(y_0, y) = \sqrt{\frac{\pi \epsilon}{a'(0)}} e^{-\frac{2}{\epsilon} A(0)}. \qquad (G.52)$$

Eq. (G.52) highlights that $I(y_0, y)$ is now independent of $y$ and $y_0$ and simply denoted $I$. Therefore Eq. (G.46) becomes

$$A(y) = K, \qquad (G.53)$$

with

$$K = -\epsilon \ln \sigma - \frac{\epsilon}{2} \ln \frac{I}{\epsilon} + \epsilon \ln x_{y_0} \qquad (G.54)$$

and the function $A(y)$ given by (G.34).

Eq. (G.53) can be expressed as a third order polynomial equation with respect to $y$ as

$$a_1 y^3 + a_2 y^2 + a_3 y + a_4 = 0 \qquad (G.55)$$

with

$$\begin{aligned} a_1 &= \frac{\zeta^2 F_1^2}{4\omega_1^2}, & a_2 &= \frac{1}{4} \left( \frac{(3\hat{\gamma}^{\mathrm{st}} - 2)\zeta^2 F_1^2}{\omega_1^2} - \alpha_1^2 \right), & a_3 &= \frac{(\hat{\gamma}^{\mathrm{st}} - 1)(3\hat{\gamma}^{\mathrm{st}} - 1)\zeta^2 F_1^2}{4\omega_1^2} - \alpha_1 K, & a_4 &= \frac{(\hat{\gamma}^{\mathrm{st}} - 1)^2\hat{\gamma}^{\mathrm{st}}\zeta^2 F_1^2}{4\omega_1^2} - K^2 \end{aligned} \qquad (G.56)$$

and $\hat{\gamma}^{\mathrm{st}}$ given by Eq. (G.7) (details are given in Appendix G.10).

Eq. (G.55) is solved analytically using Cardano's formula (see again Appendix G.10 for details). One obtains

$$\begin{aligned} \hat{y}_{\mathrm{stoch},a}^{\mathrm{dyn}} &= \begin{cases} r_0 & \text{if } \sigma < \sigma_3, \\ r_2 & \text{if } \sigma_3 < \sigma < \sigma_2 \end{cases} & (G.57a) \\ \hat{y}_{\mathrm{stoch},b}^{\mathrm{dyn}} &= \begin{cases} r_0 & \text{if } \sigma < \sigma_3, \\ r_2 & \text{if } \sigma_3 < \sigma < \sigma_2 \end{cases} & (G.57b) \end{aligned}$$

where $r_0$ and $r_2$ are given by Eq. (G.85), as depicted in Fig. G.13. The expressions of $\sigma_2$ and $\sigma_3$ are given respectively by Eqs. (G.82) to (G.84).

A second approximate expression of the stochastic dynamic bifurcation point is obtained using (G.48b) in Eq. (G.53) which becomes a second order polynomial equation with respect to $y$. Solving the latter leads to the following less accurate approximate, but easier to interpret, expression

$$\hat{y}_{\mathrm{stoch},b}^{\mathrm{dyn}} = 2 \left[ - \frac{(\hat{\gamma}^{\mathrm{st}})^{3/2} \omega_1 \epsilon}{\zeta F_1 (3 \hat{\gamma}^{\mathrm{st}} + 1)} \left( \frac{3 \ln(\hat{\gamma}^{\mathrm{st}})}{2} + \ln \left( \frac{8 \pi \omega_1}{\epsilon \zeta F_1 (3 \hat{\gamma}^{\mathrm{st}} + 1)} \right) + 4 \ln(\sigma) - 4 \ln(x_{y_0}) \right) \right]^{1/2} . \qquad (G.58)$$

Note that both in Eqs. (G.57) and (G.58) the dependence on the initial value $y_0$ is lost. However, contrary to the deterministic case, the stochastic dynamic bifurcation point depends on the slope $\epsilon$.

Deterministic and stochastic dynamic bifurcation points, given by Eqs. (G.42), (G.57) and (G.58) respectively, are plotted in Fig. G.6 as functions of the noise level $\sigma$. The deterministic bifurcation point $\hat{y}_{\mathrm{det}}^{\mathrm{dyn}}$ is plotted for two values of the initial condition $y_0$, i.e. $y_0 = -\hat{\gamma}^{\mathrm{st}}$ and $-3\hat{\gamma}^{\mathrm{st}}/4$. The stochastic bifurcation points are plotted for two values of the parameter $\epsilon$, i.e. $\epsilon = 0.002$ and $0.01$. The other parameters are given by (G.45) with in addition $x_{y_0} = 0.01$. First, we can see that the higher the noise level the closer the two approximate expressions of the stochastic dynamic bifurcation points. This is because the higher the noise level the more the bifurcation delay is reduced, hence the more the Taylor series (G.48b) is valid. Secondly, Fig. G.6 shows graphically the domains of existence of the regimes listed at the end of the Sect. G.4.1.
<figure><img src="image_133.png" /><figcaption>FIGURE G.6. Deterministic and stochastic dynamic bifurcation points, given by Eqs. (G.42), (G.57) and (G.58) respectively, as functions of the noise level $\sigma$. The deterministic bifurcation point $\hat{y}_{\text{det}}^{\text{dyn}}$ is plotted for two values of the initial condition $y_0$, i.e. $y_0 = -\hat{\gamma}^{\text{st}}$ and $-3\hat{\gamma}^{\text{st}}/4$. The stochastic bifurcation points are plotted for two values of the parameter $\epsilon$, i.e. $\epsilon = 0.002$ and 0.01. The other parameters are given by (G.45) with in addition $x_{y_0} = 0.01$.</figcaption></figure>

The following definitions are chosen for the boundary values (with respect to the noise level $\sigma$) between the different regimes previously mentioned: the boundary value between Regime I and Regime II, denoted $\sigma_{I/II}$, corresponds to the intersection between $\hat{y}_{\text{det}}^{\text{dyn}}$ and $\hat{y}_{\text{stoch,a}}^{\text{dyn}}$. This means in Fig. G.6 that if $\sigma < \sigma_{I/II}$ (Regime I), $\hat{y}_{\text{det}}^{\text{dyn}}$ is the actual dynamic bifurcation point. On the contrary, if $\sigma > \sigma_{I/II}$ (Regime II), $\hat{y}_{\text{stoch,a}}^{\text{dyn}}$ or $\hat{y}_{\text{stoch,b}}^{\text{dyn}}$ (depending on the approximation retained) is the actual dynamic bifurcation point. The boundary value between Regime II and Regime III, denoted $\sigma_{II/III}$, is the value at which $\hat{y}_{\text{stoch,a}}^{\text{dyn}}$ (and $\hat{y}_{\text{stoch,b}}^{\text{dyn}}$) no longer exists. The expressions of $\sigma_{I/II}$ and $\sigma_{II/III}$ are given in next section.

### G.4.4 Domains of existence of the regimes

The boundary value between Regime I and Regime II, with respect to the noise level $\sigma$ and denoted $\sigma_{I/II}$, is obtained solving

$$ \mathcal{D}(y) = S(y). \qquad (G.59) $$

with respect to $\sigma$. Therefore, from (G.40) and (G.52), we find

$$ \sigma_{I/II} = x_{y_0} \left( \frac{a'(0)\epsilon}{\pi} \right)^{1/4} e^{\frac{2}{\epsilon}(A(0)-A(y_0))}. \qquad (G.60) $$

The second boundary value between Regime II and Regime III, again with respect to the noise level $\sigma$ denoted $\sigma_{II/III}$, corresponds to the value of $\sigma$ for which $\hat{y}_{\text{stoch,a}}^{\text{dyn}}$ and $\hat{y}_{\text{stoch,b}}^{\text{dyn}}$ do not exist anymore<fnref n="5" /> (see Fig. G.6). The expression of $\sigma_{II/III}$ is obtained by noting that $\hat{y}_{\text{stoch,a}}^{\text{dyn}}$ and $\hat{y}_{\text{stoch,b}}^{\text{dyn}}$ vanish at $\sigma = \sigma_{II/III}$. Therefore, one may solve $\hat{y}_{\text{stoch,b}}^{\text{dyn}} = 0$ that, from (G.58), leads to

$$ \sigma_{II/III} = x_0 \left( \frac{1}{\hat{\gamma}^{\text{st}}} \right)^{3/8} \left( \frac{\epsilon \zeta F_1 (3\hat{\gamma}^{\text{st}} + 1)}{8\pi\omega_1} \right)^{1/4} \qquad (G.61) $$

with $\hat{\gamma}^{\text{st}}$ given by Eq. (G.7).

Using Eqs. (G.60) and (G.61) respectively the values of $\sigma_{I/II}$ and $\sigma_{II/III}$ corresponding to situations plotted in Fig. G.6 are given in Tab. G.2.

<footnote n="5">5. This corresponds to the fact that the solutions of Eq. (G.53) for both cubic and quadratic approximations become complex for $\sigma > \sigma_{II/III}$.</footnote>
<figure><img src="image_134.png" /><figcaption>**FIGURE G.7.** Regions of existence of the Regimes I, II and III (a) in the $(y_0, \epsilon, \sigma)$-space, (b) in the $(y_0, \sigma)$-plane for $\epsilon = 0.002$ and (c) in the plane $(\epsilon, \sigma)$-plane for $y_0 = -0.34$. The expressions of $\sigma_{I/II}$ and $\sigma_{II/III}$ are given by Eqs. (G.60) and (G.61) respectively.</figcaption></figure>

<table><caption>**TABLEAU G.2.** Values $\sigma_{I/II}$ and $\sigma_{II/III}$, computed through Eqs. (G.60) and (G.61) respectively, corresponding to situations plotted in Fig. G.6.</caption><thead><tr><th></th><th></th><th>y<sub>0</sub> = -γ<sup>st</sup></th><th>y<sub>0</sub> = -3γ<sup>st</sup>/4</th></tr></thead><tbody><tr><td rowspan="2">σ<sub>I/II</sub></td><td>ε = 0.002</td><td>1.88 · 10<sup>-14</sup></td><td>4.39 · 10<sup>-7</sup></td></tr><tr><td>ε = 0.01</td><td>1.20 · 10<sup>-5</sup></td><td>3.59 · 10<sup>-4</sup></td></tr><tr><td rowspan="2">σ<sub>II/III</sub></td><td>ε = 0.002</td><td>1.16 · 10<sup>-3</sup></td><td>1.16 · 10<sup>-3</sup></td></tr><tr><td>ε = 0.01</td><td>1.17 · 10<sup>-3</sup></td><td>1.17 · 10<sup>-3</sup></td></tr></tbody></table>

By means of Eqs. (G.60) and (G.61) one can also depict the regions of existence of each of the three regimes in the $(y_0, \epsilon, \sigma)$-space (see Fig. G.7(a)). A representation of these regions is also given in the $(y_0, \sigma)$-plane for $\epsilon = 0.002$ and in the plane $(\epsilon, \sigma)$-plane for $y_0 = -0.34$ in Figs. G.7(b) and G.7(c) respectively. One can see that the boundary between Regime II and Regime III mainly depends on $\sigma$. On the contrary, the boundary between Regime I and Regime II depends simultaneously on $\sigma$, $y_0$ and $\epsilon$. That means that the use of Eq. (G.42) (the deterministic bifurcation point) or Eq. (G.58) (the stochastic bifurcation point) as a theoretical prediction of the actual dynamic bifurcation point observed on numerical simulations does not depend only on the noise level. Fig. G.7 shows that the dependence on $\epsilon$ is the most important. Indeed, the above value $\sigma_{I/II} = 3.15 \cdot 10^{-10}$ has been obtained for $\epsilon = 0.002$. The value becomes $\sigma_{I/II} = 7.21 \cdot 10^{-17}$ for $\epsilon = 10^{-3}$ and $\sigma_{I/II} = 1.44 \cdot 10^{-69}$ for $\epsilon = 2 \cdot 10^{-4}$. In other words, if $\epsilon = 2 \cdot 10^{-4}$, to simulate the deterministic case, the number of digits of precision used by the computer must be at least equal to 69.

The latter observation shows that, in some instances, even with an extremely low level, the noise can influence the bifurcation delay. That means for example that a deterministic approach can fail to predict the bifurcation delay observed on numerical integration of the model in which the round-off errors act as an additive noise with a very low level [a6].

### G.4.5 Probability density function of the stochastic averaged amplitude

In addition to being able to obtain the expected value of the squared amplitude (G.39), Eq. (G.36) allows us to calculate the probability density function (PDF) of the amplitude $x_y$ without having to solve the Fokker-
Planck equation associated to Eq. (G.29). Indeed, a known result of stochastic calculus (see e.g. [116], Chap. 4) is that the Itô integral of a deterministic function $f(t)$, i.e. $\int_0^t f(t)dW_t$, is a Gaussian process with mean equal to zero and variance $v = \int_0^t f(t)^2dt$. Therefore, from (G.36), $x_y$ is a Gaussian process with mean $\mathbb{E}[x_y] = x^{\det}(y)$ and variance

$$ \mathbb{E}[x_y^2] - \mathbb{E}[x_y]^2 = S(y) = \frac{\sigma^2}{\epsilon} e^{\frac{2}{\epsilon} A(y)} \int_{y_0}^{y} e^{-\frac{2}{\epsilon} A(y)} dy \quad (G.62) $$

with $\int_{y_0}^{y} e^{-\frac{2}{\epsilon} A(y)} dy$ given by (G.52).

The associated Gaussian PDF is therefore

$$ \rho(x, y) = \frac{1}{\sqrt{2\pi S(y)}} e^{-\frac{(x-x_{\det}(y))^2}{2S(y)}}. \quad (G.63) $$

This means that, for a given value of $y$, the probability that $x_y$ lies between $x$ and $x + dx$ is $\rho(x, y)dx$. Then, we define the function $R(y) = \rho(x_{y_0}, y)$ for which the probability that $x_y$ lies between $x_{y_0}$ and $x_{y_0} + dx$ is $R(y)dx$. From (G.63), the expression of the function $R(y)$ is

$$ R(y) = \frac{1}{\sqrt{2\pi S(y)}} e^{-\frac{(x_{y_0} - x^{\det}(y))^2}{2S(y)}}. \quad (G.64) $$

In the Regime I we have $D(y) = x^{\det}(y)^2 \gg S(y)$ and the function $R(y)$ can be approximated by a Dirac delta function as $R(y) \approx \delta(x_{y_0} - x^{\det}(y))$. Denoting $\ell(y) = x_{y_0} - x^{\det}(y)$ the roots of which are $y_0$ and $\hat{y}_{\det}^{\text{dyn}}$ (see Eq. (G.41)) and using the composition property of Dirac delta function one obtains

$$ R(y) \approx R_I(y) = \frac{\delta(y - y_0)}{|\ell'(y_0)|} + \frac{\delta(y - \hat{y}_{\det}^{\text{dyn}})}{|\ell'(\hat{y}_{\det}^{\text{dyn}})|} \quad (G.65) $$

when the system evolves according to Regime I.

In the Regime II, because $D(y) \ll S(y)$, Eq. (G.64) becomes

$$ R(y) \approx R_{II}(y) = \frac{1}{\sqrt{2\pi S(y)}} e^{-\frac{x_{y_0}^2}{2S(y)}}. \quad (G.66) $$

We have therefore

$$ R'_{II}(y) = \frac{S'(y)(x_{y_0}^2 - S(y))}{\sqrt{2\pi S(y)^{5/2}}} e^{-\frac{x_{y_0}^2}{2S(y)}}. \quad (G.67) $$

which vanishes for $x_{y_0} = S(y)$ (it can be shown that this corresponds to a maximum of $R'_{II}(y)$) and therefore, according to the Definition G.1, for $y = \hat{y}_{\text{stoch},a}^{\text{dyn}}$.

In Fig. G.8 the function $R(y)$ is plotted using parameters (G.45) and three values of the noise level, i.e. $\sigma = 10^{-4}$, $\sigma = 10^{-6}$ and $\sigma = 10^{-11}$. For the first two values of the noise level the system is in Regime II and it is in the Regime I for the third value (see Fig. G.10 in the next section). In Figs. G.8(a) and G.8(b) one can see that the maximum of $R(y)$ is obtained for $y = \hat{y}_{\text{stoch},a}^{\text{dyn}}$ as predicted by Eq. (G.67). Eq. (G.65) is illustrated in Fig. G.8(c) in which we can see that in Regime I the function $R(y)$ becomes a Dirac delta function translated to $y = \hat{y}_{\det}^{\text{dyn}}$. In the figure only the positive values of $y$ are represented and $y_0 < 0$. Consequently, the first term of the sum in the right-hand side of Eq. (G.65) does not appear.

# G.5 COMPARISON BETWEEN THEORETICAL RESULTS AND NUMERICAL SIMULATION

## G.5.1 Comparison in term of the expected value of the squared amplitude

The aim of this section is to validate both the stochastic averaging procedure in term of expected value of the squared amplitude and the analytical predictions $\hat{y}_{\det}^{\text{dyn}}$ of the deterministic dynamic bifurcation point given by Eq. (G.42) and those of the stochastic bifurcation points, $\hat{y}_{\text{stoch},a}^{\text{dyn}}$ and $\hat{y}_{\text{stoch},b}^{\text{dyn}}$, given by Eqs. (G.57) and (G.58) respectively.

To achieve that, after getting the time series of $p_t$ from the numerical integration of (H.27) $^6$, the corresponding amplitude (denoted $A_y$) is computed as a function of $y$. To that end, the time series of $p_t$ and $y_t$ are

<footnote n="6">6. All numerical simulations of the Itô stochastic differential equations are performed using the function ItoProcess of the Wolfram Mathematica software.</footnote>
<figure><img src="image_135.png" /><figcaption>**FIGURE G.8.** The function $R(y)$ defined by Eq. (G.64). The quantities $\hat{y}_{\text{det}}^{\text{dyn}}$, $\hat{y}_{\text{stoch},a}^{\text{dyn}}$ and $\hat{y}_{\text{stoch},b}^{\text{dyn}}$ are depicted by vertical green line and vertical red dashes line respectively. The parameters (G.45) are used and (a) $\sigma = 10^{-4}$, (b) $\sigma = 10^{-6}$ and (c) $\sigma = 10^{-11}$. Moreover $y_0 = -0.34$, $x_{y0} = 0.01$ and $\epsilon = 0.002$.</figcaption></figure>

divided into $N$ intervals and on each of them the maximum of $p_t$ ($\max(p)_i$ with $i = 1, \dots, N$) and the mean of $y_t$ (mean($y_i$) with $i = 1, \dots, N$) are computed. Then $A_y$ is defined as $\max(p)_i$ as a function of mean($y_i$) for $i$ from 1 to $N$. The procedure is repeated 50 times and the expected value $\mathbb{E}[A_y^2]$ is computed over the 50 realizations. Simultaneously, the expected value $\mathbb{E}[x_y^2]$ of $x_y$ is computed over 50 realizations of the numerical integration of the averaged equation (G.29).

The results are shown in Fig. G.9 which represents : (on the left) a plot of $\sqrt{\mathbb{E}[A_y^2]}$ (dashed red line), $\sqrt{\mathbb{E}[x_y^2]}$ (blue line), the static bifurcation diagram (red) computed in Appendix G.11 and $\hat{y}_{\text{det}}^{\text{dyn}}$, $\hat{y}_{\text{stoch},a}^{\text{dyn}}$ and $\hat{y}_{\text{stoch},b}^{\text{dyn}}$ depicted by vertical black line, vertical green line and vertical magenta dashed line respectively; (on the right) same as on the left without the static bifurcation diagram, with an horizontal black dashed at $\sqrt{\mathbb{E}[x_y^2]} = x_{y0}$ and using a logarithm scale for the $\sqrt{\mathbb{E}[x_y^2]}$-axis. The set of parameters (G.45) is used and, from top to bottom : $\sigma = 10^{-4}$, $\sigma = 10^{-6}$ and $\sigma = 10^{-11}$. The initial conditions are $y_0 = -0.34$ and $x_{y0} = 0.01$.

The agreement through numerical integration between expected of the squared amplitude of the stochastic slow dynamics and the initial full order system for different situations corresponding to Regimes I and II is considered as a successful assessment of the stochastic averaging procedure.

In Fig. G.9(a) and G.9(b) the system is in Regime II. Therefore, the noise cannot be ignored and the dynamic bifurcation must be determined by $\hat{y}_{\text{stoch},a}^{\text{dyn}}$ or $\hat{y}_{\text{stoch},b}^{\text{dyn}}$. The figure shows that the latter provides a good approximation of the dynamic bifurcation point. Indeed, through the Definition G.1, the dynamic bifurcation point is the abscissa of the point of intersection between $\mathbb{E}[x_y]$ and the horizontal at $\mathbb{E}[x_y] = x_{y0}$. The logarithmic scale used in Fig. G.9 (right column) allows us to locate this intersection and to see that the vertical green line, corresponding to $\hat{y}_{\text{stoch},a}^{\text{dyn}}$, passes also through this intersection and the vertical magenta dashed line, corresponding to $\hat{y}_{\text{stoch},b}^{\text{dyn}}$, passes a little on the left (all the more as $\sigma$ is decreased). Of course the expression $\hat{y}_{\text{stoch},a}^{\text{dyn}}$ is more accurate than $\hat{y}_{\text{stoch},b}^{\text{dyn}}$ because it corresponds to a lesser degree of approximation. In
<figure><img src="image_136.png" /></figure>

**FIGURE G.9.** On the left : a plot with the expected values of $A_y^2$ (red dashed line) and $x_y^2$ (blue line) obtained using 50 realizations of the numerical integration of Eqs (H.27) and (G.29) respectively, the stable branches of the static bifurcation diagram (SBD, gray line) computed in Appendix G.11 and $\hat{y}_{\text{det}}^{\text{dyn}}$, $\hat{y}_{\text{stoch,a}}^{\text{dyn}}$ and $\hat{y}_{\text{stoch,b}}^{\text{dyn}}$ depicted by vertical black line, vertical green line and vertical magenta dashed line respectively. On the right : same as on the left without the static bifurcation diagram, with an horizontal black dashed at $\sqrt{\mathbb{E}[x_y^2]} = x_{y0}$ and using a logarithm scale for the $\sqrt{\mathbb{E}[x_y^2]}$-axis. The set of parameters (G.45) is used and, from (a) to (c), $\sigma = 10^{-4}$, $\sigma = 10^{-6}$ and $\sigma = 10^{-11}$. Moreover $y_0 = -0.34$, $x_{y0} = 0.01$ and $\epsilon = 0.002$.

Fig. G.9(c) the system is in Regime I, consequently, the noise can be neglected and the dynamic bifurcation point is determined by $\hat{y}_{\text{det}}^{\text{dyn}}$. A good agreement between theoretical prediction and numerical simulations can still be observed on the right column. Indeed, the vertical black line passes through the point of intersection between $\mathbb{E}[x_y]$ and the horizontal line at $\mathbb{E}[x_y] = x_{y0}$.

The position in the $(y_0, \sigma)$-plane of the three situations depicted in Fig. G.9 is shown in Fig. G.10 and superimposed to the regions of existence of the Regimes I, II and III.
<figure><img src="image_137.png" /><figcaption>**FIGURE G.10.** The regions of existence of the Regimes I, II and III in the $(y_0, \sigma)$-plane in which the red points show the position of the three cases depicted in Fig. G.9. The initial conditions are $y_0 = -0.34$ and $x_{y_0} = 0.01$ and three values of the noise level are used: $\sigma = 10^{-4}$, $\sigma = 10^{-6}$ and $\sigma = 10^{-11}$ which correspond to cases 1 to 3 respectively.</figcaption></figure>

### G.5.2 Comparison in term of probability density function

In this section the theoretical PDF $\rho(x, y)$ (see Eq. (G.63)) is compared to two histograms: one in built from the full order system (H.27) and the other from the stochastic slow dynamics (G.29). In both cases, 2000 realizations of the system are computed and we take, for each realization, the value of the considered random variable ($x_y$ or $p_t$) for a given value of $y$ denoted $\tilde{y}$ (chosen arbitrary to be $\tilde{y} = 0.15$). The comparison, presented in Fig. G.11, is performed for three values of the noise level (the same as chosen previously in Figs. G.8 and G.9) and shows an excellent agreement between theoretical and numerical results. In Fig. G.11(c), the system is in Regime I and, as expected, the PDF appears as a Dirac delta function. In this case the histograms are also in agreement with the theoretical PDF but this cannot be seen in the figure.

We choose a value of $\tilde{y}$ smaller than the dynamic bifurcation point but results using a larger value (not presented here) shows also a good agreement between theoretical PDF and histograms as long as the linear approximation is valid and therefore for a value $y$ not too far from the dynamic bifurcation point.

## G.6 DISCUSSION AND CONCLUSION

This article is inspired by work initially done in the applied mathematics community. The approach adopted is general and can be adapted to other systems undergoing a dynamic Hopf bifurcation.

For the clarinet model studied, the stochastic averaging procedure preserves the properties of the dynamic behavior of the original system. Thus, the same bifurcation delay is observed in the original and the averaged system. The different behaviors observed according to the amplitude of the random forcing, the initial condition or the rate of change of the bifurcation parameter are in line with those found in other fields of physics or in applied mathematics. In particular, the conclusions are the same as for another simple clarinet model already studied a few years ago. The latter model is a discrete-time model and is based on very different simplifying assumptions: it results from the discretization of a delay model which therefore retains an infinite number of acoustic modes. In contrast, the model studied in this article retains only a single acoustic mode.

Studying the dynamic bifurcations of ODE-based models is nevertheless necessary. Firstly, because these are the models used by researchers in the context of instrument making. Secondly, because recent works underline the importance of taking into account the temporal dynamics of the bifurcation parameters (in particular those controlled by the musician [41]). The study carried out in this article already allows to obtain analytically the dynamic oscillation thresholds according to the main parameters controlled by the musician. This information could be linked to important notions for the musician and the instrument manufacturer such as the ease of playing. This opens interesting perspectives such as the study of dynamic bifurcations of a more realistic model than the one studied in this article. The first improvement is undoubtedly to complete this model with other acoustic modes whose importance is already known in the static case [37].
<figure><img src="image_138.png" /><figcaption>FIGURE G.11. Comparison between theoretical PDF $\rho(x, y)$ given by Eq. (G.63) and two histograms: one built from the full order system (H.27) and the other from the stochastic slow dynamics (G.29). In both cases, 2000 realizations of the system are computed and we take, for each realization, the value of the considered random variable ($x_y$ or $p_t$) for a given value of $y$ denoted $\tilde{y}$ (chosen arbitrary to be $\tilde{y} = 0.15$). The comparison is performed for three values of the noise level: (a) $\sigma = 10^{-4}$, (b) $\sigma = 10^{-6}$ and (c) $\sigma = 10^{-11}$. Moreover $y_0 = -0.34$, $x_{y0} = 0.01$ and $\epsilon = 0.002$.</figcaption></figure>

## G.7 APPENDIX A. THE 1-DIMENSIONAL ITÔ'S FORMULA

Let the following Itô differential equation

$$dx_t = m(x_t, t)dt + \sigma(x_t, t)dW_t, \qquad (G.68)$$

where $m$ and $\sigma$ are real functions and $W_t$ is the so-called Wiener process.

Let $f(x_t, t) \in C^2(\mathbb{R}, \mathbb{R}^+)$ (i.e. $f$ is twice continuously differentiable on $(\mathbb{R}, \mathbb{R}^+)$). Then $f(x_t, t)$ is also Itô process, and

$$df(x_t, t) = \frac{\partial f}{\partial t} dt + \frac{\partial f}{\partial x} dx_t + \frac{1}{2} \frac{\partial^2 f}{\partial x^2} \sigma(x_t, t)^2 dt \qquad (G.69a)$$

$$= \left( \frac{\partial f}{\partial t} + m(x_t, t) \frac{\partial f}{\partial x} + \frac{\sigma(x_t, t)^2}{2} \frac{\partial^2 f}{\partial x^2} \right) dt + \sigma(x_t, t) \frac{\partial f}{\partial x} dW_t. \qquad (G.69b)$$

Eq. (G.69) is the 1-dimensional Itô's formula (for more details and proof see [153], Chap. 4).

## G.8 APPENDIX B. GENERAL FORMULATION OF THE STOCHASTIC AVERAGING METHOD

In this appendix the stochastic averaging method [114, 191] is briefly described. For that we consider the following system of differential equations in standard form

$$\dot{x}_t = f(x_t, t) + g(x_t, t)\eta_t \qquad (G.70)$$
where $\mathbf{x}_t \in \mathbb{R}^n$. If the deterministic vector function $\mathbf{f}(\mathbf{x}_t, t) \in \mathbb{R}^n$ and matrix function $\mathbf{g}(\mathbf{x}_t, t) \in \mathbb{R}^n \times \mathbb{R}^n$ satisfy certain requirements [114] and the elements of the vector $\eta_t$ are broadband processes, with zero means, then the slow (or averaged) dynamics of Eq. (G.70) may be approximated by the following Itô equations

$$d\mathbf{x}_t = \mathbf{m}(\mathbf{x}_t)dt + \boldsymbol{\sigma}(\mathbf{x}_t, t)d\mathbf{W}_t, \qquad (G.71)$$

where $\mathbf{W}_t \in \mathbb{R}^n$ is a vector of $n$ Wiener processes. The vector $\mathbf{m}$ and the matrix $\boldsymbol{\sigma}$ are called drift vector and diffusion matrix respectively and defined by

$$\mathbf{m} = T^{\text{av}} \left\{ \mathbf{f} + \int_{-\infty}^{0} \mathbb{E} \left[ \left( \frac{\partial(\mathbf{g}\eta)}{\partial \mathbf{x}} \right)_t (\mathbf{g}\eta)_{t+\tau} \right] d\tau \right\}, \qquad (G.72)$$

where $\frac{\partial(\mathbf{g}\eta)}{\partial \mathbf{x}}$ is the Jacobian matrix of $\mathbf{g}\eta$, and

$$\boldsymbol{\sigma}\boldsymbol{\sigma}^T = T^{\text{av}} \left\{ \int_{-\infty}^{+\infty} \mathbb{E} \left[ (\mathbf{g}\eta)_t (\mathbf{g}\eta)_{t+\tau}^T \right] d\tau \right\} \qquad (G.73)$$

where $\{\cdot\}^T$ and $\mathbb{E}\{\cdot\}$ denotes respectively the transpose and the expected value of $\{\cdot\}$. $T^{\text{av}}$ is an averaging operator defined as follows

$$T^{\text{av}} \{\cdot\} = \lim_{T \to +\infty} \frac{1}{T} \int_{t_0}^{t_0+T} \{\cdot\} dt. \qquad (G.74)$$

It should be noted that in the case of a periodic variables with period $T_0$ (which is the case in this paper), the operator $T^{\text{av}}$ becomes a classical Krylov-Bogolyubov time averaging over one period $T_0$, i.e.

$$T^{\text{av}} \{\cdot\} = \frac{1}{T_0} \int_{t_0}^{t_0+T_0} \{\cdot\} dt \qquad (G.75)$$

and the result is independent of $t_0$.

## G.9 APPENDIX C. DERIVATION OF THE EXPRESSION OF THE DETERMINISTIC DYNAMIC BIFURCATION POINT

In this appendix we give the details of the resolution of Eq. (G.41). First we state $y + \hat{\gamma}^{\text{st}} = X^2$ and then it can be shown that, using Eq. (G.34), Eq. (G.41) takes the following form

$$(X - X_0) \left( \zeta F_1 X^2 + (\zeta F_1 X_0 - \alpha_1 \omega_1) X - \zeta F_1 + \zeta F_1 X_0^2 - \alpha_1 X_0 \omega_1 \right) = 0 \qquad (G.76)$$

where $X_0^2 = y_0 + \hat{\gamma}^{\text{st}}$. Obviously $X = X_0$ and therefore $y = y_0$ is a solution of (G.76). The second term of the product in the left-hand side of (G.76) is a second order polynomial equations whose roots $X_1$ and $X_2$ are

$$X_1 = \frac{1}{2\zeta F_1} \left( \alpha_1 \omega_1 - \zeta F_1 X_0 - \sqrt{\alpha_1^2 \omega_1^2 + 2\alpha_1 \zeta F_1 X_0 \omega_1 + \zeta^2 F_1^2 (4 - 3X_0^2)} \right) \qquad (G.77)$$

and

$$X_2 = \frac{1}{2\zeta F_1} \left( \alpha_1 \omega_1 - \zeta F_1 X_0 + \sqrt{\alpha_1^2 \omega_1^2 + 2\alpha_1 \zeta F_1 X_0 \omega_1 + \zeta^2 F_1^2 (4 - 3X_0^2)} \right). \qquad (G.78)$$

The initial value $y_0$ is always chosen to be larger than $-\hat{\gamma}^{\text{st}}$ (because the mouth pressure $\gamma$ must be larger than zero). Therefore, $X$ must be larger than zero and, for a realistic set of parameters, only $X_2$ is positive. Consequently, the expression of the deterministic dynamic bifurcation point is given by

$$\hat{y}_{\text{det}}^{\text{dyn}} = X_2^2 - \hat{\gamma}^{\text{st}}. \qquad (G.79)$$
G.10 APPENDIX D. DERIVATION OF THE EXPRESSION OF THE STOCHASTIC DYNAMIC BIFURCATION POINT

Using Eq. (G.34), Eq. (G.53) takes the following explicit form

$$K = \frac{\zeta F_1 \sqrt{\hat{\gamma}^{\text{st}} + y}}{2\omega_1} (y - 1 + \hat{\gamma}^{\text{st}}) - \frac{\alpha_1 y}{2}. \qquad (G.80)$$

To obtain a solvable cubic form (i.e. without square root), Eq. (G.80) is transformed into

$$(K + \frac{\alpha_1 y}{2})^2 = \frac{\zeta^2 F_1^2 (\hat{\gamma}^{\text{st}} + y)}{4\omega_1^2} (y - 1 + \hat{\gamma}^{\text{st}})^2 \qquad (G.81)$$

which yields Eq. (G.55).

In the remaining of this appendix the Cardano's method (see e.g. [185]) is used to solve the latter, i.e. $a_1y^3 + a_2y^2 + a_3y + a_4 = 0$.

First, the following parameters are introduced

$$p = -\frac{a_2^2}{3a_1^2} + \frac{c_3}{a_1} \quad \text{and} \quad q = \frac{a_2}{27a_1} \left( \frac{2a_2^2}{a_1^2} - \frac{9a_3}{a_1} \right) + \frac{a_4}{a_1}.$$

The discriminant $\Delta$ is defined as $\Delta = -(4p^3 + 27q^2)$. Then:

1. If $\Delta < 0$, one root is real and two are complex conjugate.

2. If $\Delta = 0$, all roots are real and at least two are equal.

3. If $\Delta > 0$, all roots are real and unequal.

A typical example of the discriminant $\Delta$, plotted as a function of the noise level $\sigma$, is depicted in Fig. G.12 for a typical set of parameters. It can be shown that $\Delta$ can be expressed as a fourth order polynomial equation with respect to $\ln \sigma$ which has two distinct roots

$$\ln(\sigma_1) = \frac{2\alpha_1^3\omega_1^3 - 27\alpha_1\hat{\gamma}^{\text{st}}\zeta^2 F_1^2\omega_1}{54\zeta^2 F_1^2\omega_1\epsilon} - \frac{2(\alpha_1^2\omega_1^2 + 3\zeta^2 F_1^2)^{3/2} - 9\alpha_1\zeta^2 F_1^2\omega_1}{54\zeta^2 F_1^2\omega_1\epsilon} \\ - \frac{1}{2}\ln\left(\frac{2\sqrt{2\pi}\hat{\gamma}^{\text{st}}^{3/4}\omega_1 e^{-\frac{(\hat{\gamma}^{\text{st}}-1)\sqrt{\hat{\gamma}^{\text{st}}\zeta F_1}}{\omega_1\epsilon}}}{\sqrt{(3\hat{\gamma}^{\text{st}}+1)\zeta F_1\omega_1\epsilon}}\right) + \ln(x_{y_0}) \qquad (G.82)$$

and

$$\ln(\sigma_2) = \frac{9\alpha_1(1 - 3\hat{\gamma}^{\text{st}}) + \frac{2(\alpha_1^3\omega_1^3 + (\alpha_1^2\omega_1^2 + 3\zeta^2 F_1^2)^{3/2})}{\zeta^2 F_1^2\omega_1}}{54\epsilon} \\ - \frac{1}{2}\ln\left(\frac{2\sqrt{2\pi}\hat{\gamma}^{\text{st}}^{3/4}\omega_1 e^{-\frac{(\hat{\gamma}^{\text{st}}-1)\sqrt{\hat{\gamma}^{\text{st}}\zeta F_1}}{\omega_1\epsilon}}}{\sqrt{(3\hat{\gamma}^{\text{st}}+1)\zeta F_1\omega_1\epsilon}}\right) + \ln(x_{y_0}) \qquad (G.83)$$

and a double root

$$\ln(\sigma_3) = \frac{\alpha_1(\hat{\gamma}^{\text{st}} - 1)}{2\epsilon} - \frac{\log\left(\frac{2\sqrt{2\pi}\hat{\gamma}^{\text{st}}^{3/4}\omega_1 e^{-\frac{(\hat{\gamma}^{\text{st}}-1)\sqrt{\hat{\gamma}^{\text{st}}\zeta F_1}}{\omega_1\epsilon}}}{\sqrt{(3\hat{\gamma}^{\text{st}}+1)\zeta F_1\omega_1\epsilon}}\right)}{2} + \ln(x_{y_0}). \qquad (G.84)$$

In the example shown in Fig. G.12 one has: $\ln(\sigma_1) = -53.1$, $\ln(\sigma_2) = -6.8$ and $\ln(\sigma_3) = -26.6$.

If $\sigma_1 < \sigma < \sigma_2$ we have $\Delta > 0$ except at $\sigma = \sigma_3$ for which $\Delta = 0$. In general, when $\Delta > 0$ the three real roots are written using trigonometric functions as follows

$$r_k = 2\sqrt{\frac{-p}{3}} \cos\left(\frac{1}{3}\arccos\left(\frac{3q}{2p}\sqrt{\frac{3}{-p}}\right) + \frac{2k\pi}{3}\right) - \frac{a_2}{3a_1} \quad \text{with} \quad k = 0, 1, 2. \qquad (G.85)$$

Fig. G.13 shows the roots $r_k$ ($k = 0, 1, 2$) as functions of the noise level $\sigma$ using again the parameters (G.45) and $x_{y_0} = 0.01$. The stochastic dynamic bifurcation point, denoted $\hat{y}_{\text{stoch},a}^{\text{dyn}}$, is equal to $r_0$ if $\sigma < \sigma_3$ and to $r_2$ if $\sigma_3 < \sigma < \sigma_2$. This choice is justified by means of a numerical resolution which shows that this is the unique positive solution of Eq. (G.80).
<figure><img src="image_139.png" /><figcaption>FIGURE G.12. Discriminant $\Delta$ of (G.55) as a function the natural logarithm of the noise level $\ln \sigma$. The set of parameters (G.45) is used.</figcaption></figure>

<figure><img src="image_140.png" /><figcaption>FIGURE G.13. The roots $r_k$ ($k = 0, 1, 2$) of Eq. (G.55) and the stochastic dynamic bifurcation $\hat{y}_{\text{stoch,a}}^{\text{dyn}}$ as functions of the natural logarithm of the noise level $\sigma$. The parameters (G.45) are used and $x_{y_0} = 0.01$.</figcaption></figure>

## G.11 APPENDIX E. STATIC BIFURCATION DIAGRAM

The static bifurcation diagram is the result of the bifurcation analysis of a deterministic dynamical system with constant bifurcation parameter, here Eq. (H.26) with a constant value of $y$. It plots, as a function of the considered bifurcation parameter (here $y$), the possible steady-state regimes (fixed points and periodic motions) indicating their stability.

We use here the averaging procedure to obtain the approximated analytical bifurcation diagram of the one-mode model. For that Eq. (G.27) is considered without noise et with a constant bifurcation parameter $y$, i.e.

$$ \frac{dx_t}{dt} = F(x_t, y), \qquad (G.86) $$

where the function $F(x_t, y_t)$ is given by Eq. (G.20).

The fixed points $x^e$ of (G.86) are obtained by solving $F(x, y) = 0$. We obtain three solutions: the trivial solution $x_1^e = 0$ and two non trivial solutions, one is negative and one is positive. Only the the positive non trivial solution is retained and denoted $x_2^e$, its expression is

$$ x_2^e = 4 \sqrt{\frac{\frac{2\zeta F_1(3\hat{\gamma}^{\text{st}} + 3y - 1)}{3\omega_1\sqrt{\hat{\gamma}^{\text{st}} + y}} - \frac{4\alpha_1}{3}}{\frac{\zeta F_1(\hat{\gamma}^{\text{st}} + y + 1)}{\omega_1(\hat{\gamma}^{\text{st}} + y)^{5/2}}}} \qquad (G.87) $$
where the expression of $\hat{y}^{\text{st}}$ is given by Eq. (G.7). In the lossless case with $\alpha_1 = 0$ and $\hat{\gamma}^{\text{st}} = \frac{1}{3}$, Eq. (G.87) reduces to $x_2^e = 4\sqrt{2}(3y+1)\sqrt{\frac{y}{9y+12}}$.

The trivial solution corresponds to zero equilibrium position of (H.26) whereas the non trivial solution characterizes its periodic steady-state regimes.

As we know, the trivial fixed point $x_1^e$ is stable if $y < 0$ and unstable if $y > 0$. One can be shown that the non trivial fixed point $x_2^e$ exists only for $y > 0$ and is stable.
To be submitted to Chaos : an interdisciplinary journal of nonlinear science (2023) [p1]

## COMPUTING THE SEPARATRIX BETWEEN DYNAMIC BASINS OF ATTRACTION OF LOW DIMENSIONAL DYNAMICAL SYSTEMS WITH A TIME-VARYING PARAMETER : APPLICATION TO A MODEL OF MUSICAL INSTRUMENT

Baptiste BERGEOT, Soizic TERRIEN and Christophe VERGEZ

### ABSTRACT

We investigate the behavior of one-dimensional non-autonomous dynamical systems obtained by slowly varying over time the bifurcation parameter of the corresponding autonomous systems, i.e., whose bifurcation parameter is constant. In this context, the study focuses on the case for which the time-varying parameter crosses the bistability domain of the corresponding autonomous system. Considering the time-varying parameter as a new (slow) state variable, the considered non-autonomous one-dimensional system becomes a two-dimensional fast-slow system. In the bistability domain, the latter has attracting manifolds (resp. a repelling manifold) associated with the two stable branches (resp. unstable branch) of the bifurcation diagram of the corresponding autonomous system with constant parameter. In the framework of the geometric singular perturbation theory, we define the separatrix in the phase space between what we call the dynamic basin of attraction of the attracting manifolds. Reverse time numerical integration is used to practically compute this separatrix. Finally the proposed methodology is illustrated on a simple musical reed instrument model. This highlights which asymptotic regime is reached (played note or silence) in the bistability domain, depending on how the musician's control is varied over time.

**Keywords** : Single reed instruments, Time-varying control parameter, Global dynamic bifurcation, Bases of attraction

### LEAD PARAGRAPH

Multistability is a common phenomenon in nonlinear dynamics. It corresponds to the coexistence, in a dynamical system, of stable solutions for a given set of parameters. A multistable system can switch abruptly from a state to another when a parameter is slightly varied, a phenomenon sometimes referred to as critical transition. In this context, we investigate the behavior of one-dimensional non-autonomous dynamical systems obtained by slowly varying over time a bifurcation parameter of the corresponding autonomous systems, i.e., whose bifurcation parameter is constant. The study focuses on the case for which the time-varying parameter crosses the bistability domain of the corresponding autonomous system. Depending on the characteristics of the time-varying parameter, when the bistability domain is crossed, the system can remain close to its current position or leave it abruptly to reach another parameter dependent attractor (called attracting manifold in the paper). Considering the time-varying parameter as another (slow) state variable and then the studied model as a fast-slow system, we define and compute the separatrix in the phase space (which includes the time-varying parameter as another state variable) between what we call the dynamic basin of attraction of a given attractor, i.e., the set of initial conditions leading the system to this attractor. We illustrate our methodology using a simple bistable model of reed musical instrument.

## H.1 INTRODUCTION

Nonlinear dynamical systems can possess coexisting stable solutions (or states) for a given set of parameters. This phenomenon is called *multistability* and can imply stable states of different nature such as equilibria, periodic, quasiperiodic and chaotic solutions. To illustrate the concepts considered below we consider the following *n*-dimensional autonomous (i.e., with no explicit dependence on time) dynamical system

$$ \dot{x} = f(x, \gamma), \qquad (H.1) $$
with $x \in \mathbb{R}^n$ the vector of state variables, $\gamma$ the considered control parameter and $f$ a nonlinear function.

Multistability can be related to the phenomenon of “critical transition”, i.e., a sudden and large change in the output of the system induced by a small change in the input. The link between multistability and critical transition is discussed in the survey article [63] and applications are presented in physics, neuroscience, climate science, biology and ecology. Critical transitions are also referred to as “tipping points”, in particular in climate sciences, and can relate to different phenomena. Ashwin et al. [7] propose three categories based on the mathematical description of the underlying mechanisms : noise-induced transitions, bifurcation-induced transitions and rate-induced transitions. Noise-induced transitions, where noise make a multistable system switch from one stable state to another, are not considered in the present article. Then there are the bifurcation-induced transitions (or dynamic bifurcations [15]) and the rate-induced transitions. Both of these transitions consider $(n + 1)$-dimensional autonomous systems of the following form

$$ \dot{x} = f(x, y), \qquad (H.2a) $$

$$ \dot{y} = \epsilon g(y), \qquad (H.2b) $$

associated to Eq. (H.1) but with a time-varying bifurcation parameter $y$ which therefore becomes a new state variable. The time variation of the latter is slow compared to the original time scale of Eq. (H.1); this is accounted for here by the small dimensionless parameter $\epsilon$ in the right-hand side of Eq. (H.2b). $y$ is called the slow variable, the variables of the original state vector $x$ are called the fast variables and Eq. (H.2) is a fast-slow dynamical system. Solutions of Eq. (H.1) (both transient and stationary) are no more solutions of Eq. (H.2) since $y$ is no constant in time anymore. However, stable stationary solutions of Eq. (H.1) are associated with attracting invariant manifolds of Eq. (H.2) that evolve in its $(n + 1)$-dimensional phase space and attract orbits which are in their vicinity [62, 107]. Bifurcation-induced transitions correspond to the slow variable $y$ passing through a bifurcation point $\gamma^*$ at which the considered stable solution of Eq. (H.1) loses stability. At this point, the corresponding manifold of Eq. (H.2) ceases to be attracting and in general the transition to another attracting manifold does not appear at $\gamma^*$ but at a larger value through a phenomenon of bifurcation delay [18, 118, 198]. The rate-induced transitions, identified by Wieczorek et al. [212], are transitions between stable attracting manifolds of Eq. (H.2) without the presence of noise or a bifurcation point. In this case the critical transition is only due to the time variation of parameter under consideration. Considering equilibria of Eq. (H.1) (becoming a one-dimensional manifold of Eq. (H.2)), Ashwin et al. [7] derived the critical value of the parameter $\epsilon$ above which the orbit of the system leaves the attracting manifold it was going along.

In Ref. [6], Ashwin et al. introduce a formalism that uses so-called pullback attractors to describe the phenomena of bifurcation-induced transitions and rate-induced transitions, again considering equilibria. The method is extended by Alkhayoun and Ashwin [3] to the case of periodic solutions of Eq. (H.1), becoming two-dimensional manifolds of Eq. (H.2). In a realistic physical system, critical transitions may be associated with a combination of the three mechanisms mentioned above.

The works cited in the previous paragraph focus on the conditions causing an orbit to leave the neighborhood of a manifold that it was following. Therefore, these works deal with the beginning of the critical transition. However, the question of the end of the transition, when several other attracting manifolds can be potentially reached, is poorly addressed in the literature. In the case of dynamical system with constant parameters (i.e., of type (H.1)) for an asymptotically stable solution to be reached, the initial conditions must be in its basin of attraction, i.e., the set of initial conditions leading to it. What happens to the notion of basin of attraction when the control (or bifurcation) parameter slowly varies over time, i.e., when Eq. (H.2) is considered? In other words, what happens to an orbit of Eq. (H.2) when it leaves the neighborhood of a given attractive manifold (or a manifold which has been locally attractive) and several attractive manifolds exist?

In this paper, we propose first elements of response to these questions. To do so, (i) we consider a one-dimensional original system (i.e., $x \in \mathbb{R}$ in Eq. (H.1) and (H.2)) and (ii) we consider the framework of the geometric singular perturbation theory (GSPT) [62, 107] which provides general results on the behavior of fast-slow dynamical systems. In this context, the concepts of dynamic basin of attraction (DBA) and of the separatrix between two DBAs are introduced. Basically, these concepts can be exemplified as follows : if, in the case of a one-dimensional nonlinear ODE with a constant control parameter (i.e., of type (H.1) with $x \in \mathbb{R}$), the system is bistable with two coexisting stable equilibria whose basins of attraction, called here static basin of attraction (SBA), are separated by an unstable equilibrium; then, for the same ODE with a slowly time-varying control parameter (i.e., of type (H.2) with $x \in \mathbb{R}$) the stable equilibria becomes attracting one-dimensional manifolds that evolve in the two-dimensional phase space. Depending on the initial condition considered, one or other of these manifolds is reached by the system orbit. The set of initial conditions, in the two-dimensional phase space, leading to a given attracting manifold constitutes the DBA of this manifold. Logically, the Separatrix Between two DBA (SDBA) also becomes a one-dimensional manifold. The present research aims (i) to propose mathematical definitions of DBA and SDBA and to propose a method to compute
the SDBA and (ii) to illustrate the method on a toy model of reed instrument with a linearly time-increasing mouth pressure parameter.

Self-sustained musical instruments - such as wind instruments, bowed string instruments or singing voice - are nonlinear dynamical systems known to produce a diversity of different sound regimes. An important specificity is that sound production in a musical context is associated to a time-varying control. Indeed, control parameters are modified continuously by the instrument player. In particular, wind instruments players control air pressure in their mouth with variations over time finely tuned, first to start the sound and then to obtain the desired sound effect. However, in general, when the corresponding mathematical models of musical instruments are analyzed, the control parameters are assumed to be constant. Consequently, the instruments are modeled by autonomous nonlinear systems of differential equations (ODEs) which can have a wealth of solutions including equilibria, periodic and quasiperiodic solutions. In this context, the silence corresponds to the trivial equilibrium of the instrument model and a musical note to a periodic solution. In general, for a given set of constant control parameters, this desired periodic solution is not the only stable solution of the model. In practice, from the point of view of instrument playing, the size of the SBA of a given periodic solution can be linked to the ease with which the corresponding note can be played. For example, in the case of a reed instrument such as the saxophone, musicians know that for a given fingering corresponding to a low pitch note, the first harmonic of the desired note are sometimes easier to play. It is also on this principle of multistability that brass instruments, such as trumpets, operate.

In physical models of wind musical instruments, state variables are classically the modal pressures and their derivatives: these are the fast variables of Eq. (H.2). The considered control parameter is the pressure inside the mouth of the musician. We assume that it varies slowly over time and it is considered as another state variable: the slow variable such as $y$ for Eq. (H.2). In this context, previous theoretical works be the authors [a5, a6, a15] investigated the emergence of oscillations in a simple models of reed instruments in which the pressure in the mouth, which the control parameter, varies over time. This allowed to interpret experimental results [a4] as the manifestation of the bifurcation delay phenomenon. After studying this in the case of the appearance of oscillations, we investigate here the influence of the temporal dynamics of the mouth pressure on the nature of the regime reached in the presence of multistability. A problem that has been preliminary studied by Colinot et al. [43] through a posteriori observations on results of direct time numerical integration of a saxophone model. Authors showed that different regimes may be reached by modifying the mouth pressure dynamics. Here, we focused on the underlying phenomenons allowing us to interpret these kinds of observations. Although the dynamics observed in this particular system corresponds to a bifurcation-induced transition, the concepts of BDA and SDBA might also de of interest in the case of rate-induced transitions.

The paper is organized as follows. In Section H.2, the single reed instrument model under consideration is briefly presented. The static bifurcation diagram of the model is computed in Section H.3. Section H.4 provides elements of geometric singular perturbation theory needed for the mathematical definition of the DBA and of the separatrix between two of them. The latter are defined in Section H.5.1 in which the method used to compute the separatrix is also presented. The method is then illustrated on the studied model in Section H.5.2. A discussion about the case of a blowing pressure which saturates in the bistability domain and about the extension of the method to more realistic models of instruments is proposed in Section H.6. Finally, concluding remarks and some perspectives are given in Section H.7.

## H.2 BACKGROUND ON THE PHYSICAL MODEL OF REED INSTRUMENT

The model of single reed instrument model considered in the paper consists in a classical toy model obtained under strong assumptions (see Section H.8.1). We focus here on dynamic basins of attraction mentioned above and describe a method to define and compute the separatrix between two basins. For sake of conciseness, details on the clarinet model (which is classical in musical acoustics) are presented in the Appendix H.8.

The state variables are the pressure $p$ at the entrance to the instrument resonator (i.e., inside the mouthpiece) and its first time derivative and the bifurcation parameter is the blowing pressure $\gamma$. We consider the averaged dynamics of this model, in which the amplitude $x$ of $p$ is a state variable, with a time-varying blowing pressure (the control parameter under consideration in this work), now considered as another state variable and denoted $y$ (see Section H.8.2).

Equations of the averaged dynamics are of the form of Eq. (H.2), i.e.,

$$ \dot{x} = f(x, y), \qquad (H.3a) $$

$$ \dot{y} = \epsilon g(y), \qquad (H.3b) $$

with $(x, y) \in \mathbb{R}^2$. Eq. (H.3a) is the equation of the physical instrument model and Eq. (H.3b) describes the
**TABLEAU H.1.** Modal and embouchure parameters. Taken from Ref. [81]. Unless stated otherwise, these parameter values are used through out the article.

<table><tr><td>Modal parameters</td><td>ε<sub>1</sub> = 1/36.6<br/>ω<sub>1</sub> = 1440<br/>Z<sub>1</sub> = 50</td></tr><tr><td>Embouchure parameter</td><td>ζ = 0.1</td></tr></table>

<figure><img src="image_141.png" /><figcaption>**FIGURE H.1.** Static bifurcation diagram of the clarinet obtained by averaging procedure. A bistability domain is observed between $\gamma = 1.11$ and $\gamma = 1.34$ and depicted by a blue colored area. The parameters given in Table H.1 are used.</figcaption></figure>

time variation of the blowing pressure $y$. The function $f$ is derived from the so-called nonlinear characteristic (NLC) of the instrument exciter which links the volume flow through the reed channel to the pressure difference between the mouth and the mouthpiece and has been widely studied in the literature [50, 98].

In the following sections, we refer to *static* when the mouth pressure is time-constant (i.e. $\dot{y} = 0$) and to *dynamic* when the mouth pressure is time-varying (i.e. $\dot{y} = \epsilon g(y)$).

## H.3 STATIC BIFURCATION DIAGRAM OF THE MODEL

The *static bifurcation diagram* is a commonly used tool in musical acoustics (see e.g. Refs. [37, 81]). It provides useful information beyond the purely static case, and has shown its interest in the dynamic case, in particular in a model of reed musical instrument [a5, a6, a15]. Here, the static bifurcation diagram is computed by stating $\dot{y} = 0$ with $y = \text{const.} = \gamma$ in Eq. (H.3). Fixed points of $\dot{x} = f(x, \gamma)$ (with $f$ given by Eq. (H.36)) are computed by solving $f(x, \gamma) = 0$. In practice, we use the function NSolve of the Mathematica software [216].

The trivial fixed point is the equilibrium solution of the original (non-averaged) dynamics (H.27) and the non trivial fixed points correspond to the periodic solution of (H.27). The stability of the fixed points is determined subsequently from the sign of $\partial_x f(x, \gamma)$.

Fig. H.1 shows the static bifurcation diagram obtained for the modal parameters of the air-column and the embouchure parameter $\zeta$ given in Table H.1. Unless stated otherwise, the parameter values given in Table H.1 are used through out the article. Fig. H.1 shows that the model is bistable between $\gamma = 1.11$ and $\gamma = 1.34$. Within this range of the parameter value the trivial fixed point and a non trivial fixed point are stable and, lying between them, the other non trivial fixed point is unstable.

For a given stable solution, the *static basin of attraction* (SBA) is the set of initial conditions leading to this solution. Here, for our 1D system $\dot{x} = f(x, \gamma)$, the separatrix between the SBAs of the two stable fixed points in the bistability domain is trivial, and is the unstable fixed point.

## H.4 ELEMENTS OF GEOMETRIC SINGULAR PERTURBATION THEORY

In this section some basics of the geometric singular perturbation theory (GSPT) [107] are recalled. System (H.3) is considered for sake of illustration. In the GSPT framework (i) the latter is a (1, 1)-fast-slow system, i.e., a dynamical system with one fast variable ($x$, here the amplitude of the mouthpiece pressure) and
<figure><img src="image_142.png" /><figcaption>FIGURE H.2. Times series $x(\tau)$ (green) and $p(\tau)$ (blue) obtained from the numerical integration of the averaged system (H.4) and of the original system (H.27) (replacing $\gamma$ by $y(t')$ and switching from $t'$ to $t$, see Appendix H.8), respectively. Both are obtained with $g(y) = 1$. The time profile of the linearly increased blowing pressure $y(\tau)$ is also represented (red). The parameters given in Table H.1 are used with $\epsilon = 0.0183$. The initial conditions are: $x(0) = p(0) = 0.1$, $\dot{p}(0) = 0$ and $y(0) = 0.15$.</figcaption></figure>

one slow variable (here $y$, the blowing pressure) and (ii) $t$ represents the fast timescale. Most of these elements of GSPT are reproduced from Berglund and Gentz [22, 24] and Kuehn [40].

Eq. (H.3) is now written on the slow timescale $\tau = \epsilon t$, Eq. (H.3) becomes

$$ \epsilon x' = f(x, y) \qquad (H.4a) $$

$$ y' = g(y), \qquad (H.4b) $$

where $\{\}' = d_{\tau}\{\}$ (the first time derivative with respect to the slow time $\tau$).

As an illustration, Fig. H.2 shows an example of possible time series $x(\tau)$ and $p(\tau)$ obtained from the numerical integration of the averaged system (H.4) and of the original system (H.27) (replacing $\gamma$ by $y(t')$ and switching from $t'$ to $t$, see Appendix H.8), respectively. We considered a linear time variation of the slow variable, i.e.,

$$ y(\tau) = \tau + y_0, \qquad (H.5) $$

or $y(t) = \epsilon t + y_0$, with $y_0 = y(0)$. Therefore, the function $g(y)$ is simply

$$ g(y) = 1. \qquad (H.6) $$

We used the parameters given in Table H.1.

The time profile of the linearly increasing blowing pressure $y(\tau)$ is also represented. Overall, the figure shows an excellent agreement between time series $x(\tau)$ and $p(\tau)$ which validates the averaging procedure.

The time evolution of slow-fast systems is characterized by possible successions of fast epochs and slow epochs. This is shown in Fig. H.2. Indeed, the variable $x$ first decreases rapidly to zero (first fast epoch). Then during the increase of the mouth pressure, $x$ follows zero (first slow epoch) and increases rapidly at $\tau \approx 0.65$ (second fast epoch). A slow evolution follows corresponding to oscillations of the mouthpiece pressure $p$ (second slow epoch). Finally, $x$ decreases again rapidly to zero at $\tau \approx 1.45$ (third fast epoch) and follows zero slowly (third slow epoch).

In the GSPT framework, the so-called slow subsystem is obtained by stating $\epsilon = 0$ in Eq. (H.4), which leads to

$$ 0 = f(x, y) \qquad (H.7a) $$

$$ y' = g(y). \qquad (H.7b) $$

The dynamics of (H.7) approximates the dynamics of (H.3) (or (H.4)) during slow epochs. It is restricted to the critical manifold $M_0$ defined below.
**DEFINITION H.1.** The critical manifold is defined as

$$ \mathcal{M}_0 = \{(x, y) \in \mathbb{R}^2 | x = x^*(y)\} \qquad (H.8) $$

with $x^*(y)$ the branches of the solution of

$$ f(x^*(y), y) = 0 \qquad (H.9) $$

in each interval of $y$ where $\partial_x f(x, y)$ does not vanish and therefore where we can write that $f(x, y) = 0$ is equivalent to $x = x^*(y)$ by the *implicit function theorem*.

Note that the critical manifold does not depend on the function $g(y)$.

Points on the critical manifold are equilibria (or fixed points) of the so-called *fast subsystem* defined by

$$ \dot{x} = f(x, y) \qquad (H.10a) $$

$$ \dot{y} = 0 \qquad (H.10b) $$

which describes the dynamics of the variable $x$ for a constant bifurcation parameter. This corresponds to the static case considered in Section H.3.

**DEFINITION H.2.** Let $a^*(y)$ be the linearization of the fast vector field (H.10) on $\mathcal{M}_0$, i.e., at $x = x^*(y)$, defined as

$$ a^*(y) = \partial_x f(x^*(y), y) \qquad (H.11) $$

A value $x^*(y)$ of the fast variable $x$ is a hyperbolic equilibrium point of (H.10) if $a^*(y) \neq 0$; stable if $a^*(y) < 0$ and unstable if $a^*(y) > 0$. Then, the critical manifold is called *attracting* (resp. *repelling*) if $a^*(y) < 0$ (resp. $a^*(y) > 0$) for $y \in I$ with $I$ a subset of $\mathbb{R}$. A subset $\mathcal{M}_0^{\text{nh}}$ of the critical manifold $\mathcal{M}_0$ is *normally hyperbolic* if for each point $P = (x^*(y), y) \in \mathcal{M}_0^{\text{nh}}$ we have $a^*(y) \neq 0$, i.e., $\mathcal{M}_0^{\text{nh}}$ is either attracting or repelling$^a$.

A simple form of the Fenichel's theorem (from [24]) is given below. It states that all orbits starting near an attracting branch of the critical manifold actually converge to an invariant manifold.

**THEOREM H.1 (FENICHEL 1979).** If the critical manifold $\mathcal{M}_0$ is normally hyperbolic (i.e. attracting or repelling), then there exists a manifold $\mathcal{M}_\epsilon$, which is $\mathcal{O}(\epsilon)$-close$^a$ to $\mathcal{M}_0$ and invariant under the flow (H.3)$^b$ (or (H.4)). The manifold $\mathcal{M}_\epsilon$ is normally hyperbolic, i.e., it attracts or repels neighboring orbits exponentially fast in directions normal to itself.

For a more complete statement of the Fenichel's theorem see e.g. Chap. 3 of Ref. [40].

**DEFINITION H.3.** The manifold $\mathcal{M}_\epsilon$, as obtained by the Theorem H.1, is called the *slow manifold*.

Here and as often, we say “the” manifold $\mathcal{M}_\epsilon$. Formally, this is incorrect. Indeed, the complete statement of the Fenichel's theorem stipulates that $\mathcal{M}_\epsilon$ is usually not unique. However, it also stipulates that all these manifolds are exponentially close to each other, i.e., they lie at a distance$^1 \mathcal{O}(e^{-K/\epsilon})$ from each other for some $K > 0$ with $K = \mathcal{O}(1)$. Consequently, we adopt the convention to refer to $\mathcal{M}_\epsilon$ as “the” slow manifold of Eq. (H.4).

Based on these elements of GSPT, the definition of what we referred to as a *dynamic basin of attraction* is formalized in Section H.5.

<footnote>a. 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 text below of Ref. [40]).</footnote>

<footnote>a. To be understood in the sense of the Hausdorff distance (see e.g. Chap. 3 of Ref. [40]).</footnote>

<footnote>b. $\mathcal{M}_\epsilon$ is called invariant (locally in $I$) under the flow, if $(x_0, y_0) \in \mathcal{M}_\epsilon$ implies that $(x(t), y(t)) \in \mathcal{M}_\epsilon$ as long as $y(t) \in I$ holds.</footnote>

<footnote n="1">1. Again, to be understood in the sense of the Hausdorff distance.</footnote>
## H.5 SEPARATRIX BETWEEN DYNAMIC BASINS OF ATTRACTION

This section presents the main results of this work. We give in Section H.5.1 definitions for general (1, 1)-fast-slow systems. The results are then illustrated in Section H.5.2.

### H.5.1 Definitions and method for computing separatrix between DBAs

The definitions stated in this section are valuable for (1, 1)-fast-slow systems in the form of Eq. (H.3) whose critical manifold $M_0$ (see Definition H.1) has a bistability domain. The definition of the latter is given below.

**DEFINITION H.4 (BISTABILITY DOMAIN OF $M_0$).** An open subset

$$D =]y_1, y_u[ \qquad (H.12)$$

of $\mathbb{R}$ is called a *bistability domain* of the critical manifold if for $y \in D$ Eq. (H.9) has two attracting branches $x_1^*(y)$ and $x_2^*(y)$, i.e., $a_1^*(y) = \partial_x f(x_1^*(y), y) < 0$ and $a_2^*(y) = \partial_x f(x_2^*(y), y) < 0$. Necessarily, Eq. (H.9) has also a repelling branch $x_3^*(y)$ for which $a_3^*(y) = \partial_x f(x_3^*(y), y) > 0$. The critical manifold is non normally hyperbolic at $y = y_1$ and $y = y_u$.

The bistability domain of the critical manifold $M_0$ is of course, in the case of a (1, 1)-fast-slow system, identical to the bistability domain of the corresponding static problem (see Fig. H.1), one has therefore $y_1 = 1.11$ and $y_u = 1.34$.

We can now define what we call a *dynamic basin of attraction*.

**DEFINITION H.5 (DYNAMIC BASIN OF ATTRACTION).** In the bistability domain $D$ (see Definition H.4) the critical manifold $M_0$ of (H.4) has two coexisting attracting branches $M_{0,a_1}$ and $M_{0,a_2}$ defined as

$$M_{0,a_i} = \{x \in \mathbb{R}, y \in D \mid x = x_i^*(y)\}, \quad i = 1, 2 \qquad (H.13)$$

separated by a repelling branch $M_{0,r}$

$$M_{0,r} = \{x \in \mathbb{R}, y \in D \mid x = x_3^*(y)\} \qquad (H.14)$$

Through Fenichel's Theorem H.1, Eq. (H.4) has two attracting invariant manifolds $M_{\epsilon,a_i}$ ($\epsilon$)-close of $M_{0,a_i}$ ($i = 1, 2$) and a repelling slow invariant manifold $M_{\epsilon,r}$ ($\epsilon$)-close of $M_{0,r}$. Then the *dynamic basin of attraction* (DBA) of a given attracting branch $M_{\epsilon,a_i}$ ($i = 1, 2$) is the subset of the phase space for which orbits originating from initial conditions in the DBA end up following $M_{\epsilon,a_i}$ ($i = 1, 2$) when the bistability domain $D$ is crossed.

We state in Definition H.5 “end up following” because some transient regimes can exist before the branches $M_{\epsilon,a_i}$ ($i = 1, 2$) are reached.

The separatrix between two DBAs and its nature are finally defined.

**DEFINITION H.6 (SEPARATRIX BETWEEN DBAS).** The separatrix between DBAs (SDBA) is the boundary in the phase space that separates two DBAs.

**DEFINITION H.7 (NATURE OF THE SDBA).** The SDBA is the repelling invariant manifold $M_{\epsilon,r}$ considered beyond the lower bound $y_1$ of the bistability domain until a point of the phase space having the smallest physically relevant value of the slow variable $y$.

The previous definitions are especially relevant when the increasing blowing pressure does not saturate in the bistability domain. In such a case, the system has no fixed points (at least it has none for $y < y_u$) and we want to know when the bistability zone is crossed whether or not the original (i.e., non averaged) system oscillates significantly. In the case of musical instrument, this corresponds to the fact that a sound is produced (or not) during the transient. The case of a blowing pressure that saturates in the bistability domain is discussed in Section H.6.1.

In view of the definitions H.6 and H.7, we are interested in computing the repelling invariant manifold $M_{\epsilon,r}$. Through the Fenichel's theorem H.1, if $x^*(y)$ is a normally hyperbolic branch of the critical manifold (either
attracting or repelling), the invariant manifold $M_\epsilon$ (either attracting or repelling) admits a parametric equation of the form $x = \bar{x}(y, \epsilon)$, where $\bar{x}(y, \epsilon) = x^*(y) + O(\epsilon)$. The function $\bar{x}(y, \epsilon)$ is the solution of $\epsilon \frac{dx}{dy} = \frac{f(x,y)}{g(x,y)}$ obtained by dividing (H.4a) by (H.4b). Asymptotic analytical expressions of $\bar{x}(y, \epsilon)$, in the form of a power series in $\epsilon$, can be easily obtained in regions in which the critical manifold is hyperbolic (see e.g. Ref. [24]). These approximations fail at non hyperbolic points, making them unusable for our study $^2$. Indeed, an SDBA is the extension of a repelling invariant manifold $M_{\epsilon,r}$ through non hyperbolic points of the critical manifold (see Definition H.6). The curve $\bar{x}(y, \epsilon)$ can be defined as a particular solution of Eq. (H.3) or Eq. (H.4). In the manner of Berglund and Landon [26], this particular solution, i.e., the SDBA, is here approximated numerically by means of a time reversal procedure, i.e., by a numerical integration of

$$ \epsilon x' = -f(x, y) \qquad (H.15a) $$

$$ y' = -g(x, y) \qquad (H.15b) $$

choosing an initial condition $(x_0, y_0)$ on the upper boundary of the bistability domain $D$ and on the repelling part $M_{0,r}$ of the critical manifold, i.e.,

$$ (x_0, y_0) = (x_3^*(y_u), y_u). \qquad (H.16) $$

We use the fact that, through Theorem H.1, $M_{0,r}$ is $O(\epsilon)$-close to the corresponding repelling invariant manifold $M_{\epsilon,r}$ which is attracting for (H.15). The numerical integration is then stopped at the required value of the slow variable $y$.

The time reversal method has also the advantage to be easily extended to higher dimensional systems (see Section H.6.2).

## H.5.2 Illustration of SDBAs in the averaged clarinet model with slowly time varying blowing pressure

The attracting and repelling branches of the critical manifold $M_0$ of Eq. (H.4) with the parameters given in Table H.1 are respectively the stable and unstable branches of the bifurcation diagram depicted in Fig. H.1 (see Definition H.1 and the beginning of Section H.3). The bistability domain of $M_0$, as described by Definition H.4, is $D = [y_1 = 1.11, y_u = 1.34]$. In the latter, the DBAs of each attracting invariant manifold $M_{\epsilon,a_1}$ and $M_{\epsilon,a_2}$ (see Definition H.5) are denoted DBA$_1$ and DBA$_2$, respectively. Without loss of generality, we chose arbitrarily that $M_{\epsilon,a_1}$ is the branch corresponding to $x_1^*(y) = 0$. The separatrix between DBA$_1$ and DBA$_2$ is obtained through Definition H.7. It is the extended slow repelling invariant manifold $M_{\epsilon,r}$ (associated to the repelling branch $M_{0,r}$ of the critical manifold) beyond the lower bound $y_1$ of the bistability domain to a point of the phase space having the smallest possible value of the slow variable $y$. In system (H.15) with a linear time variation of $y$, i.e., $g(y) = 1$, $x(\tau)$ diverges as $\tau$ tends to infinity. Consequently, the SDBA is approximated by an orbit of Eq. (H.15) having the appropriate initial condition given by (H.16). Here that yields $(x_0, y_0) = (x_3^*(y_u) = 0.95, y_u = 1.34)$. As mentioned in Section H.5.1, this orbit is obtained by a numerical simulation. The simulation is stopped when an arbitrary large value of $x(\tau)$ is reached. The resulting times series are shown in Fig. H.3(a) (red solid and dashed lines for $x(\tau)$ and $y(\tau)$, respectively).

Fig. H.3(b) shows the time integration of the direct time system (H.4) with an initial condition (IC) first in DBA$_1$ (in blue) and then in DBA$_2$ (in green). Both of these initial conditions are chosen close to the SDBA. In the case of $IC \in DBA_1$, we observe that $x(\tau)$ first rapidly approaches zero and remains equal to zero over time. In the case of $IC \in DBA_2$, similarly to the example in Fig. H.2, after a first similar start transient, $x(\tau)$ remains close to zero for a moment and then moves away from it, taking a non-zero value for a period of time before falling back to zero. The first (resp. second) situation seems to correspond well to the case where the trajectory tends to follow $M_{\epsilon,a_1}$ (resp. $M_{\epsilon,a_2}$) when the bistability domain $D$ is crossed. This appears more clearly in Fig. H.3(c) in which the orbits previously computed are represented in the phase plane and superimposed to the critical manifold. The stream plot of the vector field Eq. (H.4) is also shown. The latter shows, for $y < y_u$, that the orbits starting with an initial condition in the DBA$_1$ (resp. DBA$_2$), i.e., below (resp. above) the SDBA, end up following $M_{\epsilon,a_1}$ (resp. $M_{\epsilon,a_2}$).

To better understand what happens in the vicinity of the $y$-axis, a logarithmic scale is used in Fig. H.3(d) for the ordinates. This example shows that very small amplitudes are encountered. In a real-world problems, residual noise would prevent reaching such low amplitudes. In the future, stochastic studies (similar to what was done in [a15]) would be relevant to complete this work.

<footnote n="2">2. The behavior of a (1, 1)-fast-slow system near a non hyperbolic point of the critical manifold can be analytically studied by means of the methodology developed by Berglund et al. [17, 22] which consists in determining the scaling law of the system, i.e., the dependence of the system in the parameter $\epsilon$ in the neighborhood of the non hyperbolic point. This is not the purpose of the present paper because this analytical method cannot be extended for higher dimensional phase spaces.</footnote>
<figure><img src="image_143.png" /><figcaption>**FIGURE H.3.** Illustration of DBA and SDBA for a linear time-varying blowing pressure. (a) Times series of $x(\tau)$ (solid line) and $y(\tau)$ (dashed line) obtained by numerical integration of the reverse time system (H.15) with $g(y) = 1$ and an initial condition (IC) given by (H.16) (red). (b) Time series of the direct time system (H.4) with an IC first in DBA$_1$ (in blue) and in DBA$_2$ (in green). (c) The corresponding orbits (same colors are used) in the phase plane and superimposed to the critical manifold. The orbit corresponding to the reverse time integration corresponds to the SDBA. The stream plot of the vector field Eq. (H.4) with streamlines colored in light blue is also shown. (d) The same as in (c) but with a logarithmic scale for the ordinates. The parameters are given in Table H.1 with an additional $\epsilon = 0.0183$.</figcaption></figure>

The transient regimes mentioned after Definition H.5 are visible in the stream plot in Fig. H.3(c). These transient regimes are all the longer as the initial conditions are chosen close to the SDBA. This begins the following question: will these transients be short enough for the system to reach the branches $M_{\epsilon,a_i}$ before the end of the bistability domain? This question relates to the link between the duration of transients and the distance between the initial conditions and the SDBA. This will be explored in future works.

Examples similar to those of Figs. H.3(c) and H.3(d) are shown in Fig. H.4 for two other values of the parameter $\epsilon$: $\epsilon = 0.00366 < 0.0183$ (see Fig. H.4(a)) and $\epsilon = 0.0732 > 0.0183$ (see Fig. H.4(b)). This shows that the smaller $\epsilon$, the shorter (in terms of $y$) the above-mentioned transients, and the more visible the distinction between slow and fast phases. In Fig. H.4(b) in particular (top, blue line) the orbit is still relatively far from the $M_{\epsilon,a_2}$ (the $y$-axis) when the upper bound of the bistability domain is reached. This highlights the limits of the approach for large values of $\epsilon$. Nevertheless, the approach allows to interpret these borderline situations.
<figure><img src="image_144.png" /><figcaption>FIGURE H.4. The same figures as Figs H.3(c) and H.3(d) but with (a) $\hat{\epsilon} = 10^{-4}$ and therefore $\epsilon = 0.00366$ and (b) $\hat{\epsilon} = 2 \cdot 10^{-3}$ and therefore $\epsilon = 0.0732$.</figcaption></figure>

# H.6 DISCUSSION

## H.6.1 Case of a blowing pressure saturating in the bistability domain

In this section the case where the blowing pressure tends towards a finite limit $y_{targ}$ in the bistability domain is explored. More precisely, an exponential growth of the blowing pressure is considered using $g(y) = y_{targ} - y$ (with $y_{targ} \in D$) which yields

$$y(\tau) = y_{targ} + (y_0 - y_{targ})e^{-\tau}, \qquad (H.17)$$

with again $y_0 = y(0)$. Written in the slow timescale $\tau$, the system takes the following form

$$\epsilon x' = f(x, y) \qquad (H.18a)$$

$$y' = y_{targ} - y. \qquad (H.18b)$$

This differs from the case studied so far, in which the pressure does not saturate. Here fixed points are solutions of $f(x, y) = 0$ and $y_{targ} - y = 0$, that leads to $f(x, y_{targ}) = 0$ and $x = x_i^*(y_{targ})$ ($i = 1, 2, 3$). Therefore, the points $P_i^* = (x_i^*(y_{targ}), y_{targ})$ ($i = 1, 2, 3$) are the fixed points of Eq. (H.18). The Jacobian matrices $J(P_i^*)$ of (H.18) evaluated at each fixed point $P_i^*$ are

$$J(P_i^*) = \begin{pmatrix} \frac{a_i^*(y_{targ})}{\epsilon} & \frac{\partial_y f(x_i^*(y_{targ}), y_{targ})}{-1} \\ 0 & \epsilon \end{pmatrix} \qquad (H.19)$$

whose eigenvalues are $\mu_1 = \frac{a_1^*(y_{targ})}{\epsilon}$ and $\mu_2 = -1$. Because $a_1^*(y_{targ}) < 0$, $a_2^*(y_{targ}) < 0$ and $a_3^*(y_{targ}) > 0$ (see Definition H.4), $P_1^*$ and $P_2^*$ are nodes (i.e., $\mu_1\mu_2 > 0$) and $P_3^*$ is a saddle (i.e., $\mu_1\mu_2 < 0$). These classical definitions of nodes and saddle are taken from Ref. [179].

The problem therefore becomes a classic problem of finding the basins attraction of stable fixed points of a 2-dimensional system. In this case, the DBAs of the attractive branches $M_{\epsilon,a_i}$ ($i = 1, 2$) and the SBAs of
the fixed points $P_1^*$ and $P_2^*$ are identical. The separatrix between these basins is the stable manifold of the unstable fixed point $P_3^*$ [38]. The latter is again computed using a time reversal simulation. The following equation

$$ \epsilon x' = -f(x, y) \qquad (H.20a) $$

$$ y' = -y_{\text{targ}} - y. \qquad (H.20b) $$

is numerically integrated from initial condition chosen as a small perturbation of $P_3^*$.

An example is shown in Fig. H.5. Figures H.5(a) and H.5(b) show similar representations as in Figs H.3(c) and H.3(d), respectively. In Fig. H.5(a) the orange vertical line indicates the target value $y_{\text{targ}}$ of the mouth pressure. Intersections between this vertical line and the critical manifold correspond to fixed points of (H.18). Importantly, the critical manifold is independent of the function $g(y)$ and therefore identical to the one previously obtained in Section H.5.2 for a linear growth of $y$. Still in Fig. H.5(a), a zoom near the unstable fixed point $P_3^*$ is performed to highlight its stable manifold interpreted here as a SDBA (represented by a red line in the figure). We can see in Fig. H.5(b) that the SDBA intersects the x-axis for a value of $y$ just below 0.01.

Two orbits obtained by numerical integration of the direct time system Eq. (H.18) with (i) a IC first in DBA$_1$ (in blue) and (ii) in DBA$_2$ (in green) also appear in Fig. H.5. The corresponding time series are plotted in Figs. H.5(c) and H.5(d), respectively. In the case of $IC \in \text{DBA}_1$, the orbit stays close to the y-axis for a moment and then moves away from it, follows the branch $M_{\epsilon,a_1}$ of the critical manifold before stopping at $P_1^*$. For $IC \in \text{DBA}_2$, it is the point $P_2^*$ which ends up being reached. Here too transients are observed before the fixed points are reached.

## H.6.2 DBA and SDBA to more realistic models

This work introduces the notions of DBA and SDBA for (1, 1)-fast-slow systems. However, more realistic models of reed instruments cannot be reduced to (1, 1)-fast-slow systems. Working in the plane makes it easier to illustrate and understand these concepts. However, the concepts of DBA and SDBA can be generalized to $(m, n)$-fast-slow systems. Indeed, the SDBA is a particular solution of the considered system. However, there are two things to bear in mind when the fast subsystem is more than one-dimensional. Firstly, the critical manifold is also more than one-dimensional and it can be attracting (if all eigenvalues of the Jacobian matrix of the fast subsystem evaluated on the critical manifold have negative real parts), repelling (if all the eigenvalues have positive real parts) and saddle-type (if there are eigenvalues with negative real parts and others with positive real parts). Previous definitions are taken from Ref. [40] (Chap. 3). In this case of a saddle-type invariant manifold, the SDBA cannot be obtained by reversing the time, as presented above for a purely repelling manifold. Indeed, a saddle-type manifold in direct time remains saddle-type in reverse time. In this case methods of continuation of orbits could be used, for example using the software AUTO [57]. The second thing to consider is that the fast subsystem can have other types of than equilibria; such as periodic or quasiperiodic solutions. The complete bifurcation diagram of the fast subsystem (including periodic solutions and if possible quasiperiodic solutions) is not anymore the critical manifold of the studied $(m, n)$-fast-slow system which includes only the equilibrium solutions. In a $(m, n)$-fast-slow system, invariant manifolds associated to quasiperiodic solutions of the fast subsystem are more complicated to deal with, both analytically and computationally. However, periodic solutions can be treated in the same way as equilibrium solutions. First, because theoretically Berglund [22] proposes an analogue of Fenichel's theorem on the existence of an invariant manifold tracking families of periodic orbits and also because periodic solutions are easy to compute numerically. As in the case of equilibrium solutions, for the time reversal method to be used, the invariant manifold considered for the SDBA calculation must be repelling. For periodic solutions, that means that $m - 1^3$ of the Floquet multipliers associated with the periodic solution under consideration have modulus strictly larger than one (see e.g. Ref. [179], Chap. 7, for details about local stability of periodic solutions).

We provide here an illustration for the original system (H.27) (replacing $\gamma$ par $y(t')$ and switching from $t'$ to $t$, see Appendix H.8), a (2, 1)-fast-slow system (globally 3-dimensional) since the 1-dimensional complex Eq. (H.33) is a 2-dimensional real equation.

As mentioned above, the trivial fixed point of the fast subsystem (H.10) associated to the averaged dynamics (H.3) is the equilibrium solution of the original non-averaged dynamics (H.27) and non trivial fixed points correspond to periodic solutions of (H.27). In the bistability domain the fast subsystem associated to the original dynamics thus has a stable equilibrium solution, a stable and an unstable periodic solutions. The manifolds corresponding to these periodic solutions are represented in Fig. H.6 by gray and black surfaces in

<footnote n="3">3. A periodic solution of an $m$-dimensional fast subsystem is associated to $m$ Floquet multipliers but one is always unity (see Lemma 7.3 of Ref. [179]).</footnote>
<figure><img src="image_145.png" /><figcaption>**FIGURE H.5.** Illustration of DBA and SDBA for an exponential time-varying blowing pressure. (a) Orbit obtained by the numerical integration of the reverse time system (H.15) with $g(y) = y_{\text{targ}} - y$ and an initial condition (IC) chosen as a small perturbation of the non trivial fixed point $P_3^*$ (in red). This represents the SDBA of Eq. (H.18). Two orbits obtained by the numerical integration of direct time system Eq. (H.18) : with a IC first in DBA$_1$ (in blue) and in DBA$_2$ (in green). The corresponding time series are plotted in Figs. H.5(c) and H.5(d), respectively. The stream plot of the vector field Eq. (H.4) with streamlines colored in light blue is also shown is also shown. (b) The same as in (a) but with a logarithmic scale for the ordinates. The parameters are given in Table H.1 with in addition $\epsilon = 0.0732$, $y_{\text{targ}} = 1.25$.</figcaption></figure>

the $(y, p, \dot{p})$-space, respectively. In practice these manifolds, which also represent the static bifurcation diagram of the original dynamical system, are deduced from the critical manifold of Eq. (H.3). Indeed, in the $(p, \dot{p})$-plane periodic solutions of the non-averaged dynamics (H.27) is a circle whose radius is equal to the value of the corresponding fixed point of Eq. (H.3). The manifolds of the periodic solutions of the latter are then obtained by rotation of the critical manifold of Eq. (H.3) around the y-axis.

The SDBA is here the extension of a 2-dimensional repelling slow invariant manifold associated to the unstable periodic solution of the fast subsystem. In a 2-dimensional system an unstable periodic solution is associated to two Floquet multipliers, one is unity$^{??}$ and the other has modulus larger than one. Therefore, the associated invariant manifold is repelling and the reverse time procedure can be used.

In Fig. H.6 the SDBA is shown in red in the $(y, p, \dot{p})$-space. Although the SDBA is 2-dimensional, here we computed only one solution of the time-reversed system for initial conditions on the periodic solution at $y = y_u$. Here, the period of oscillations is short enough for the two-dimensional SDBA to be represented correctly by a single periodic orbit.

Results of time-domain simulation of the (2,1)-fast-slow original system (in direct time) with two different sets of initial conditions, one in DBA$_1$ (in blue) and the second one in DBA$_2$ (in green) are also shown in Fig. H.6. The parameters and IC are the same as the same as in Fig. H.3 (same parameters and initial conditions) and similar observations can be made : when $\text{IC} \in \text{DBA}_1$, the orbit first winds around the y-axis and remains in the vicinity of it. Conversely, when $\text{IC} \in \text{DBA}_2$, after a first similar start transient, the orbits stays on vicinity of the y-axis for a while and then moves away to reach the manifold associated to the stable
<figure><img src="image_146.png" /></figure>

**FIGURE H.6.** Illustration of DBA and SDBA for the original system (H.27) with a linear time-varying blowing pressure. The manifolds associated to the stable and unstable periodic solutions of the fast subsystem of the original non-averaged dynamics (H.27) are depicted by gray and black surfaces, respectively. Orbit obtained by the numerical integration of the reverse time system associated to Eq. (H.27) (replacing $\gamma$ by $y(t')$ and switching from $t'$ to $t$, see Appendix H.8) and an IC on the circle representing the periodic solution in the $(p, \dot{p})$-plane at $y = y_u$ (in red). This represents the SDBA of the 3-dimensional original dynamics. Two orbits obtained by the numerical integration of direct time system : with an IC first in DBA₁ (in blue) and in DBA₂ (in green). Same parameters as in Fig. H.3.

periodic solution. The simulation is carried out over a shorter period of time than in Fig. H.3. This is why we do not observe the orbit falling back to zero.

## H.7 CONCLUSION

In this paper the behavior of a class of one-dimensional non-autonomous dynamical systems has been investigated. These systems are obtained by slowly varying over time the bifurcation parameter of the corresponding autonomous systems (called fast subsystem in the paper) whose bifurcation parameter is constant. More precisely, the case considered corresponds to the time-varying parameter crossing a bistability domain of the associated autonomous system. The proposed methodology has been illustrated on a simple musical reed instrument model.

In the bistability domain, the non-autonomous system has attracting manifolds (resp. a repelling manifold) associated with the two stable branches (resp. unstable branch) of the bifurcation diagram of the corresponding autonomous system. The concept of dynamic basin of attraction (DBA) of a given attracting manifold has been introduced and defined as the set of initial conditions in the phase space (in which the time-varying parameter is considered as another state variable) from which the orbits end up following this attracting manifold. The separatrix between two DBAs (SDBA), i.e., the boundary in the phase space that separates the two DBAs, has been defined as the repelling invariant manifold mentioned above and considered beyond the lower bound of the bistability domain until a point of the phase space having the smallest possible value of the time-varying bifurcation parameter. In practice, the separatrix has been computed using a numerical reverse time integration with relevant initial conditions. As a first perpective, keep working on low dimension systems, it would be interesting to develop analytical methods that would make it possible to relate the characteristics of the time-varying parameter to the nature of the observed regime. Moreover, in the purpose of real-life applications, the method should be extended to more complex dynamic systems. First, the notions of DBA and SDBA remain valid for high dimensional systems. However, if the solutions (stable and unstable) of the multistable autonomous system are more complex than equilibria, the challenge will be to detect the unstable solutions to then deduce the possible SDBAs of the non-autonomous system. Secondly, multistable systems are known to be very sensitive to noise. Therefore, the influence of noise on SDBAs should be studied in the future. Finally, we believe that the general nature of the method suggests that the results obtained go beyond
musical acoustics. Application in physics, neuroscience or climate science can be envisaged, in particular to deal with the so-called critical transition problems.

# H.8 APPENDIX A. DERIVATION OF THE MODEL

## H.8.1 Classical single reed instrument model

Sound production by single reed instruments is classically modeled through the nonlinear coupling between two linear components [14, 37, 66] : the reed and the air-column inside the instrument. While blowing air through the reed channel into the instrument, the instrumentalist provides a quasi-static source of energy. The instrument and the player constitute an autonomous dynamical system. When the trivial equilibrium solution of this system becomes unstable, a sound is produced [67, 181, 215].

Since the lowest resonance frequency of the reed is one order of magnitude higher than the sound frequency for many notes, the reed is often modeled as a lossless stiffness spring [10, 155]. Therefore its position relative to rest is proportional to the pressure drop across the reed, i.e., the pressure difference between the mouth and the mouthpiece of the instrument. The linear pressure response of the air column P to the volume flow U through the reed channel is given in the frequency domain through the input impedance of the air column Z

$$P(\omega) = Z(\omega)U(\omega), \qquad (H.21)$$

where $\omega$ is the angular frequency. We point out that the model presented in the appendix already considers a dimensionless pressure and flow (see Chap. 9 of Ref. [37] for more details on the model).

The contribution at the input of the instrument of the (infinite) series of modes of the air column is taken into account in $Z(\omega)$. For computational reasons, the series is truncated to $N$ modes, where $N$ is an integer :

$$Z(\omega) = \sum_{n=1}^{N} Z_n \frac{j \epsilon_n \omega_n \omega}{\omega_n^2 + j \epsilon_n \omega_n \omega - \omega^2}, \qquad (H.22)$$

with $Z_n$, $\omega_n$ and $\epsilon_n$ the modal parameters, respectively the modal factor, the resonance angular frequency and the inverse of the quality factor of the $n^{th}$ peak of the impedance (corresponding to the $n^{th}$ mode of the air column). Eq. (H.22) can be written in the time domain

$$d_{t''t''}p_n + \epsilon_n \omega_n d_{t''}p_n + \omega_n^2 p_n = Z_n \epsilon_n \omega_n d_{t''}u, \quad \forall n \in [1, N], \qquad (H.23)$$

in which $t''$ is used to denote the original timescale, $u$ is the inverse Fourier transform of $U$ and $p_n$ is such that $p = \sum_{i=1}^{N} p_n$, where $p$ is the inverse Fourier transform of $P$ [181] and corresponds to the time evolution of the mouthpiece pressure.

Through Bernoulli's principle, the volume flow through the reed channel $u$ is related nonlinearly to the reed channel opening and the pressure difference between the mouth and the mouthpiece [50, 98]

$$u = \hat{F}(p, \gamma) = \zeta(1 + p - \gamma) \sqrt{|\gamma - p|} \operatorname{sgn}(\gamma - p) H(1 + p - \gamma) \qquad (H.24)$$

where $H$ is the Heaviside function, $\gamma$ is the dimensionless pressure in the mouth of the musician and $\zeta$ a dimensionless parameter accounting for many embouchure parameters. The parameters $\gamma$ and $\zeta$ are the control (or bifurcation) parameters of the model. The relation (H.24) is called the nonlinear characteristic (NLC) of the instrument exciter. In this work a polynomial fitting of the function $\hat{F}(p, \gamma)$ (see Eq. (H.24)), denoted $F(p, \gamma)$, is obtained using the function InterpolatingPolynomial of the Mathematica software [216] (see Fig. H.7).

Using the function $F(p, \gamma)$, Eq. (H.23) can be written using only the pressure $p$ as follows

$$d_{t''t''}p_n + \epsilon_n \omega_n d_{t''}p_n + \omega_n^2 p_n - Z_n \epsilon_n \omega_n d_{t''}p_n \partial_p F(p, \gamma) = 0, \quad \forall n \in [1, N]. \qquad (H.25)$$

A minimal model of a reed instrument including a single mode of the air-column is obtained by stating $N = 1$. In this case (H.25) becomes

$$d_{t''t''}p + \epsilon_1 \omega_1 d_{t''}p + \omega_1^2 p - Z_1 \epsilon_1 \omega_1 d_{t''}p \partial_p F(p, \gamma) = 0 \qquad (H.26)$$

In this case ($N = 1$), since $p_1 = p$, $p_1$ is replaced by $p$ in Eq. (H.26). This is clearly a minimal yet useful model of sound production in reed instruments. Indeed, it takes into account the two main control parameters, $\gamma$ and $\zeta$, adjusted by the musician and describes the physical mechanism through which sound emerges from the
<figure><img src="image_147.png" /><figcaption>FIGURE H.7. Comparison between the function $\hat{F}(p, \gamma)$ given by Eq. (H.24) (blue line) and its polynomial fitting $F(p, \gamma)$ (green line) as functions of the pressure difference $\gamma - p$.</figcaption></figure>

trivial equilibrium (i.e. silence) when a resonance of the air column is excited by an incoming flow. We finally introduce the dimensionless time $t' = \omega_1 t''$, so Eq. (H.26) takes the form of the following self-excited oscillator

$$d_{t't'}p + \epsilon_1 h(p, d_{t'}p, \gamma) + p = 0. \qquad (H.27)$$

where

$$h(p, d_{t'}p, \gamma) = d_{t'}p(1 - Z_1 \partial_p F(p, \gamma)). \qquad (H.28)$$

## H.8.2 Averaged dynamics with a slowly time-varying control parameter

Eq. (H.27), in the phase space $(p, d_{t'}p)$, is transformed to a slowly varying system using the following complex representation introducing a new variable $\xi$ as

$$\xi e^{jt'} = d_{t'}p + jp \qquad (H.29)$$

with $j^2 = -1$. Combining Eq. (H.29) and its complex conjugate yields the expressions of $p$ and $d_{t'}p$ as functions of $\xi$

$$p = \frac{\xi e^{jt'} - \xi^* e^{-jt'}}{2j}, \qquad (H.30)$$

$$d_{t'}p = \frac{\xi e^{jt'} + \xi^* e^{-jt'}}{2}, \qquad (H.31)$$

where $\xi^*$ is the complex conjugate of $\xi$. Deriving Eq. (H.29) with respect to $t'$ and using Eq. (H.31) give $d_{t't'}p$ as

$$d_{t't'}p = d_{t'}\xi e^{jt'} + j\xi e^{jt'} - \frac{j}{2}(\xi e^{jt'} + \xi^* e^{-jt'}). \qquad (H.32)$$

Substituting Eqs. (H.30), (H.31) and (H.32) into Eq. (H.27) yields the following complex non-autonomous system

$$d_{t'}\xi = -\epsilon_1 h \left( \frac{\xi e^{jt'} - \xi^* e^{-jt'}}{2j}, \frac{\xi e^{jt'} + \xi^* e^{-jt'}}{2}, \gamma \right) e^{-jt'}. \qquad (H.33)$$

Then, stating $\xi = xe^{j\varphi}$ and $\varphi = t' + \varphi$ and separating real and imaginary parts of each side of Eq. (H.33) gives

$$d_{t'}x = -\epsilon_1 h(x \sin \varphi, x \cos \varphi, \gamma) \cos \varphi, \qquad (H.34a)$$

$$d_{t'}\varphi = 1 + \epsilon_1 h(x \sin \varphi, x \cos \varphi, \gamma) \frac{\sin \varphi}{x}. \qquad (H.34b)$$

Since $0 < \epsilon_1 \ll 1$, Eq. (H.34) is simplified by means of an averaging method (see Chap. 7 of Ref. [172]) to

$$d_{t'}x = \epsilon_1 f(x, \gamma), \qquad (H.35a)$$
$$d_{t'}\phi = 1 \qquad (H.35b)$$

where $x$ and $\phi$ are uncoupled and

$$f(x, \gamma) = -\frac{1}{2\pi} \int_{0}^{2\pi} h(x \sin \alpha, x \cos \alpha, \gamma) \cos \alpha d\alpha. \qquad (H.36)$$

Details on the averaging procedure are provided in Section H.9.

The integral in Eq. (H.36) is generally hard to solve analytically, except if the function $h$ is assumed to have a polynomial form. This is the case here using the polynomial fitting $F(p, \gamma)$ of the NLC.

The dynamics of the mouth pressure (now denoted $y$) is finally added as $d_{t'}y = \hat{\epsilon}g(y)$. Then, assuming $\hat{\epsilon}/\epsilon_1 \ll 1$, switching the time from $t'$ to $t = \epsilon_1 t'$ and using the notation $\{\cdot\} = d_t\{\cdot\}$, the following system is obtained

$$\dot{x} = f(x, y) \qquad (H.37a)$$

$$\dot{y} = \epsilon g(y) \qquad (H.37b)$$

where $\epsilon$ is defined as

$$\epsilon = \frac{\hat{\epsilon}}{\epsilon_1} \qquad (H.38)$$

which is the relevant small parameter to use for this model.

## H.9 APPENDIX B. DETAILS ON THE AVERAGING PROCEDURE

Omitting here, for the sake of clarity, the dependence in $\gamma$, and introducing

$$X(x, \phi) = -h(x \sin \phi, x \cos \phi, \gamma) \cos \phi, \qquad (H.39a)$$

$$\Omega(x, \phi) = h(x \sin \phi, x \cos \phi, \gamma) \frac{\sin \phi}{x}, \qquad (H.39b)$$

Eq. (H.34) becomes

$$d_{t'}x = \epsilon_1 X(x, \phi), \qquad (H.40a)$$

$$d_{t'}\phi = 1 + \epsilon_1 \Omega(x, \phi). \qquad (H.40b)$$

The latter is a (1, 1)-fast-slow systems with $t'$ the fast timescale, $x$ the slow variable and $\phi$ the fast variable. The principle of the average procedure (see Chap. 7 of Ref. [172]), consists in assuming the following form for the variable $x$

$$x = v + \epsilon_1 u(v, \phi) \qquad (H.41)$$

with $u(v, \phi)$ a function $2\pi$-periodic in $\phi$.

Using Eq. (H.41) and

$$d_{t'}u = (1 + \epsilon_1 \Omega(v + \epsilon_1 u, \phi)) \partial_{\phi} u + d_{t'} v \partial_{v} u, \qquad (H.42)$$

Eq. (H.40a) becomes

$$d_{t'}v(1 + \epsilon_1 \partial_v u) + \epsilon_1 (1 + \epsilon_1 \Omega(v + \epsilon_1 u, \phi)) \partial_{\phi} u = \epsilon_1 X(v + \epsilon_1 u, \phi). \qquad (H.43)$$

Each $\epsilon_1$-dependent function in (H.43) is expanded in a first-order Taylor series around $\epsilon_1 = 0$, that yields

$$d_{t'}v + \epsilon_1 (1 + \epsilon_1 \Omega(v, \phi) + \mathcal{O}(\epsilon_1^2)) (1 - \epsilon_1 \partial_v u + \mathcal{O}(\epsilon_1^2)) \partial_{\phi} u = (\epsilon_1 X(v, \phi) + \mathcal{O}(\epsilon_1^2)) (1 - \epsilon_1 \partial_v u + \mathcal{O}(\epsilon_1^2)) \qquad (H.44)$$

which simplifies into

$$d_{t'}v = \epsilon_1 (X(v, \phi) - \partial_{\phi} u) + \mathcal{O}(\epsilon_1^2) \qquad (H.45)$$

Then one assumes a function $u(v, \phi)$ which satisfies the following equation

$$\partial_{\phi} u = X(v, \phi) - f(v). \qquad (H.46)$$

Integrating (H.46) with respect to $\phi$ from 0 to $2\pi$ and because $u(v, \phi)$ is $2\pi$-periodic in $\phi$ we obtain

$$f(v) = \frac{1}{2\pi} \int_{0}^{2\pi} X(v, \alpha) d\alpha, \qquad (H.47)$$
and then $f(v)$ is the average of $X(v, \phi)$ as defined in Eq. (H.36).

Substituting (H.46) into (H.45) and adding Eq. (H.40b) we obtain a simplified dynamical system with the following form

$$d_{t'}v = \epsilon_1 f(v) + \mathcal{O}(\epsilon_1^2) \qquad (H.48a)$$

$$d_{t'}\phi = 1 + \mathcal{O}(\epsilon_1). \qquad (H.48b)$$

Finally, consistently with $0 < \epsilon_1 \ll 1$, we have $\epsilon_1 f(v) + \mathcal{O}(\epsilon_1^2) \approx \epsilon_1 f(v)$, $1 + \mathcal{O}(\epsilon_1) \approx 1$ and $v \approx x$ and therefore Eq. (H.35) is proven.
TROISIÈME PARTIE

NOTICE PERSONNELLE
<figure><img src="image_148.png" /></figure>
## CURRICULUM VITÆ

### Baptiste BERGEOT

né le 14 janvier 1985 à Rennes (35)

Adresse électronique : baptiste.bergeot@insa-cvl.fr

Site internet académique : https://baptistebergeotacademic.github.io/AcademicWebsite/

**Position actuelle :** Maître de conférences (section CNU 60), INSA Centre Val de Loire<fnref n="1" />, Laboratoire de Mécanique Gabriel Lamé (LaMé, EA 7494), Équipe Dynamique interactions vibrations Structures (DivS).

## 5.1 CURSUS ACADÉMIQUE

**2015** Qualification Maître de Conférences, Section 60.

**2013** Doctorat d'acoustique, Université du Mans.

* Préparé au Laboratoire d'Acoustique de l'Université du Mans (LAUM, UMR 6613).

* Soutenue le 10 octobre 2013.

* Encadrement : Bruno Gazengel (PU, Univ. Le Mans), Christophe Vergez (DR CNRS, LMA<fnref n="2" />) et André Almeida (MCF, Univ. Le Mans).

* Intitulé : Naissance des oscillations dans les instruments de type clarinette à paramètre de contrôle variable.

* Lien TEL : https://tel.archives-ouvertes.fr/tel-01149513.

**2010** Master d'Ingénierie Mécanique et Acoustique, parcours Acoustique, Université du Mans.

**2007** Maîtrise de Sciences Physiques, parcours Physique, Université de Rennes 1.

**2003** Baccalauréat Scientifique.

## 5.2 PARCOURS PROFESSIONNEL

**Depuis 2015** Maître de Conférences

INSA Centre Val de Loire

**2013 → 2015**
**Post-doctorant**
Encadrement : Bruno Cochelin (PU, ECM) et Sergio Bellizzi (DR CNRS, LMA)

École Centrale Marseille (ECM), LMA

<footnote n="1">1. INSA Centre Val de Loire, 3 rue de la chocolaterie, CS 23410, 41034 Blois cedex, FRANCE.</footnote>

<footnote n="2">2. LMA ≡ Laboratoire de Mécanique et d'Acoustique (UMR 7031, Marseille, FRANCE).</footnote>
# 5.3 ENCADREMENTS DOCTORAUX ET SCIENTIFIQUES

## 5.3.1 Doctorats

<table><thead><tr><th>Doctorant</th><th>Informations</th></tr></thead><tbody><tr><td>Israa Zogheib</td><td>Intitulé : Comportement dynamique de systèmes mécaniques lents-rapides stochastiques<br>Date de début : 1er novembre 2023 (thèse en cours)<br>Taux de participation : 50%<br>Autre encadrant : Nils Berglund (directeur de thèse, taux de participation 50%)</td></tr><tr><td>Chérif Snoun</td><td>Intitulé : Contrôle passif des vibrations des systèmes mécaniques à l'aide d'absorbeurs dynamiques non linéaires avec prise en compte d'incertitudes paramétriques<br>Date de début : 1er octobre 2016<br>Date de soutenance : 16 juillet 2020<br>Taux de participation : 50%<br>Autre encadrant : Sébastien Berger (directeur de thèse, taux de participation 50%)</td></tr><tr><td>Duc Thinh Kieu</td><td>Intitulé : Méthodes de réduction de modèles pour l'analyse du comportement dynamique des structures complexes à paramètres incertains<br>Date de début : 1er septembre 2015<br>Date de soutenance : 10 février 2020<br>Taux de participation : 25%<br>Autres encadrants : Sébastien Berger (directeur, taux de participation 25%), Jean-Mathieu Mencik (co-directeur, taux de participation 25%) et Marie-Laure Gobert (co-encadrante, taux de participation 25%)</td></tr></tbody></table>

## 5.3.2 Stages de recherche de Master

<table><tbody><tr><td>2023</td><td>Stagiaire : Helmi Chaabene<br>Intitulé : Atténuation passive de résonance vibratoire à l'aide d'un absorbeur non linéaire bistable : étude numérique<br>Taux de participation : 100%</td></tr><tr><td>2022</td><td>Stagiaire : Joris Nguyen<br>Intitulé : Atténuation des cycles limites d'un modèle non linéaire d'éolienne à l'aide d'absorbeurs non linéaires de vibrations<br>Taux de participation : 100%</td></tr><tr><td>2021</td><td>Stagiaire : Riccardo Gega<br>Intitulé : Contrôle passif des vibrations indésirables dans les éoliennes à l'aide d'absorbeurs dynamiques non linéaires<br>Taux de participation : 100%</td></tr><tr><td>2020</td><td>Stagiaire : Joulia Salloum<br>Intitulé : Calculs stochastiques pour le dimensionnement d'un échangeur de chaleur<br>Taux de participation : 50%<br>Autre encadrant : Quentin Serra (INSA CVL, taux de participation 50%)</td></tr><tr><td>2019</td><td>Stagiaire : Barthélemy Jeandel</td></tr></tbody></table>
<table><tr><td></td><td>Intitulé : Vers la conception d'une sourdine pour instrument à anche basée sur le concept de pompage énergétique<br/>Taux de participation : 50%<br/>Autre encadrant : Christophe Vergez (LMA, taux de participation 50%)</td></tr><tr><td>2019</td><td>Stagiaire : Corentin Coustham<br/>Intitulé : Calculs stochastiques pour le dimensionnement d'un échangeur de chaleur<br/>Taux de participation : 50%<br/>Autre encadrant : Quentin Serra (INSA CVL, taux de participation 50%)</td></tr><tr><td>2016</td><td>Stagiaire : Ayda Halouani<br/>Intitulé : Calculs mécaniques robustes pour le dimensionnement de réservoirs d'hydro-gène<br/>Taux de participation : 50%<br/>Autre encadrant : Benoît Magnain (INSA CVL, taux de participation 50%)</td></tr></table>

## 5.4 LISTE DES CONTRATS DE RECHERCHE

### 5.4.1 Contrats de recherche évalués à la suite d'appels à projets

* **[Sure-Methonum]** (Sureté et ingénierie système : Méthode numériques pour les approches probabilistes) : Projet régional LAVOISIER (2015-2016).

Partenaires : INSA CVL, CEA.

Montant : 22k€.

Rôle : participant au contrat par l'encadrement d'un stage de Master 2.

* **[DSG DPT2MA]** (Déploiement de Plateformes Technologiques en Matériaux, Mécanique et Acoustique ultrasonore) : Projet de Recherche et d'innovation dans le cadre du Dialogue Stratégique et de Gestion 2021-2022 de l'INSA CVL.

Partenaires : INSA CVL, LaMé, GREMAN.

Montant : 135k€.

Durée : 2 ans.

Rôle : co-responsable scientifique pour le LaMé et encadrant d'un stage de Master 2.

**Descriptif du projet.** Collaboration avec le laboratoire GREMAN afin d'associer au contrôle passif des vibrations à l'aide d'absorbeurs dynamiques non linéaires, la récupération d'énergie à partir d'un dispositif piézoélectrique. Cette collaboration s'inscrit dans un projet de recherche et d'innovation lauréat dans le cadre du Dialogue Stratégique et de Gestion 2021-2022 de l'INSA CVL. L'objectif est d'étudier la capacité d'un absorbeur non linéaire couplé à un dispositif piézoélectrique de récupération d'énergie à réduire les vibrations néfastes à une structure primaire tout en récupérant cette énergie. L'idée est d'augmenter le transfert d'énergie de la structure mécanique primaire à protéger vers le matériau piézoélectrique, à l'aide de l'absorbeur non linéaire. Un démonstrateur expérimental est en cours de conception et de réalisation dans le cadre d'un stage de Master 2.

* **[BQR DiSPaREVMT]** (Dispositifs Semi-Passif de Récupération d'Énergie Vibratoire dans les Machines Tournantes réceptrices de type éolienne basé sur le concept de pompage énergétique et la piézoélectricité) : Bonus Qualité Recherche (BQR) 2023-2024 de l'INSA Centre Val de Loire.

Partenaires : INSA CVL, LaMé, GREMAN.

Montant : 15k€.

Durée : 1 an.

Rôle : porteur et responsable scientifique pour le LaMé.
**Descriptif du projet.** Le travail de recherche imaginé pour ce projet s'inscrit dans la continuité du projet DSG [DSG DPT2MA] décrit plus haut. Le but serait d'adapter les résultats obtenus dans ce dernier aux problèmes d'instabilités vibratoires dans les machines tournantes réceptrices de type éolienne. Le problème sera abordé à la fois d'un point de vue théorique/numérique (développement d'une méthode de conception du NES) et expérimental (mise en place d'expériences de laboratoire reproduisant les instabilités de flottement ou de couplage rotor/mat puis leur contrôle à l'aide d'un NES).

* **[COREVE]** (Contrôle et Récupération d'Énergie Vibratoire dans les éoliennes)

Il s'agit d'un appel à « projets de recherche d'initiative académique » (APR IA) de la région Centre-Val de Loire remporté en juillet 2023 et qui débutera début 2024.

**Titre long** : Dispositifs passifs pour le contrôle et la récupération de l'énergie vibratoire nuisible dans les systèmes mécaniques tournants : application aux éoliennes

Partenaires : INSA CVL, LaMé, GREMAN.

Montant : 113 k€.

Durée : 3 ans.

Rôle : responsable scientifique

**Descriptif du projet.** (cf. Section 4.2.1.1).

### 5.4.2 Contrats de recherche de gré à gré

* **[Caline 1]** (Calcul incertains d'échangeurs de chaleur 1) : Contrat industriel de recherche (2019).

Partenaires : Air Liquide, LaMé.

Montant : 14k€.

Durée : 1 an.

Rôle : participant au contrat par l'encadrement d'un stage de Master 2.

* **[Caline 2]** (Calcul incertains d'échangeurs de chaleur 2) : Contrat industriel de recherche (2020).

Partenaires : Air Liquide, LaMé.

Montant : 14k€.

Durée : 1 an.

Rôle : participant au contrat par l'encadrement d'un stage de Master 2.

## 5.5 RESPONSABILITÉS

### 5.5.1 Responsabilités scientifiques

* Membre du Conseil de Laboratoire du Laboratoire de Mécanique Gabriel Lamé (LaMé, EA 7494)

Période d'exercice : 20 février 2020 - 31 décembre 2023

* Référent sur le campus de Blois de l'INSA Centre Val de Loire du Master *mention Mécanique, Génie Civil, Matériaux, Structures* adossé au Laboratoire de Mécanique Gabriel Lamé. Le Master est co-accrédité par 3 établissements : l'Université de Tours, l'Université d'Orléans et l'INSA Centre Val de Loire

Période d'exercice : depuis la rentrée universitaire 2018 (rentrée de la première promotion du Master)
### 5.5.2 Responsabilités pédagogiques

* Responsable de l'option *Ingénierie Mécanique et Conception*, option de 5ᵉ année du département *Génie des Systèmes Industriels* de l'INSA Centre Val de Loire.

Période d'exercice : depuis la rentrée universitaire 2018

## 5.6 EXPERTISE SCIENTIFIQUE

### 5.6.1 Relecture d'articles

Relecteur de 32 articles pour les journaux suivants : Aerospace Science and Technology (3 article), Automatica (1 article), CEAS Aeronautical Journal (1 article), Chaos, Solitons & Fractals (1 article), International Journal of Non-Linear Mechanics (4 articles), Journal of Computational and Nonlinear Dynamics (1 article), Journal of Engineering Mechanics (1 article), Journal of Sound and Vibration (5 article), Journal of Vibration and Acoustics (1 article), Journal of Vibration and Control (1 article), Mechanism and Machine Theory (1 article), Nonlinear Dynamics (10 articles), The Journal of the Acoustical Society of America (1 article), The Journal of the Acoustical Society of America Express Letters (1 article).

### 5.6.2 Participation à des jurys de thèse

* Examinateur de la thèse de Rémy Tanays intitulée « *Contrôle passif non linéaire des instabilités dynamiques* », préparée à l'ISAE<fnref n="3" /> et soutenue le 12 décembre 2022 à Toulouse.

## 5.7 RAYONNEMENT ET DIFFUSION SCIENTIFIQUE

### Séminaires et communications dans des journées thématiques

1. Journées annuelles du GdR ExModeli<fnref n="4" />, FEMTO-ST Besançon

Titre : *Atténuation passive de vibrations auto-entretenues au moyen d'un absorbeur dynamique non linéaire bistable*

Date : 10 novembre 2023

2. Séminaire de l'Institut Denis Poisson (IDP UMR CNRS 7013, Orléans, France)

Titre : *Systèmes lents-rapides en mécanique vibratoire : application au contrôle passif non linéaire de vibrations auto-entretenues et à l'étude des phénomènes transitoires dans les instruments de musique à anche*

Date : 23 mars 2023

3. Séminaire du Laboratoire de Mécanique et d'Acoustique (LMA UMR CNRS 7031, Marseille, France)

Titre : *Atténuation passive d'oscillations auto-entretenues à l'aide d'absorbeurs dynamiques non linéaires de type NES : approches analytiques et numériques*

Date : 7 février 2023

4. Séminaire du Laboratoire de Mécanique des Structures et des Systèmes Couplés (LMSSC, EA 3196, CNAM, Paris, France)

Titre : *Atténuation passive d'oscillations auto-entretenues à l'aide d'absorbeurs dynamiques non linéaires de type NES : approches analytiques et numériques*

Date : 12 octobre 2022

<footnote n="3">3. ISAE ≡ Institut Supérieur de l'Aéronautique et de l'Espace</footnote>

<footnote n="4">4. Groupe de Recherche « Exploitation Et Modélisation Des Dynamiques Non Linéaires ».</footnote>
5. Séminaire du Laboratoire de Mécanique et d'Acoustique (LMA UMR CNRS 7031,
Marseille, France)

Titre : *Atténuation d'instabilités vibratoires à l'aide d'absorbeurs dynamiques non linéaires de type « NES » : prédiction théorique des régimes d'oscillations*

Date : 15 janvier 2019

6. Journées annuelles du GdR DYNOLIN<fnref n="5" />, CNAM Paris

Titre : *Steady-state regimes of a helicopter ground resonance model including a nonlinear energy sink*

Date : 14 octobre 2014

7. Séminaire du Laboratoire de Mécanique et d'Acoustique (LMA UMR CNRS 7031,
Marseille, France)

Titre : *Naissance des oscillations dans les instruments de type clarinette à paramètre de contrôle variable*

Date : 12 novembre 2013

<footnote n="5">5. Groupe de Recherche « Dynamique non linéaire ».</footnote>
# ACTIVITÉS D'ENSEIGNEMENT

## 6.1 INTRODUCTION

J'ai eu la chance de participer à la formation de jeunes scientifiques et ingénieurs depuis mon doctorat à l'Université du Mans. Mes activités d'enseignement se situent dans les domaines de la mécanique, de l'acoustique et des vibrations au niveaux Licence et Master. Elles sont décrites ci-dessous, classées par établissement où elles ont été dispensées : à l'INSA Centre Val de Loire et à l'Université de Tours en tant que Maître de Conférences et à l'Université du Mans en tant que doctorant.

Certains des supports (notes de cours et/ou diaporamas) des cours mentionnés ci-dessous sont disponibles au téléchargement sur mon site internet académique :

https://baptistebergeotacademic.github.io/AcademicWebsite/SWBB_teaching.html

## 6.2 À L'INSA CENTRE VAL DE LOIRE

### Vibrations non linéaires (2019-aujourd'hui)

Ce cours est une création. Il comporte 6h de Cours Magistraux et 6h de Travaux Dirigés pour les étudiants en 5e année FISE<fnref n="1" /> du département GSI<fnref n="2" /> option IMC<fnref n="3" /> et les étudiants du Master Mécanique, Génie Civil, Matériaux, Structures adossé au Laboratoire de Mécanique Gabriel Lamé.

L'objectif du cours est de donner aux étudiants une introduction aux méthodes d'analyses des vibrations non linéaires. Le Chapitre 1 permet de mettre en contexte le cours en présentant d'abord quelques généralités sur les non-linéarités en mécanique. Il présente ensuite, à l'aide d'exemples, certaines propriétés propres aux systèmes vibratoires non linéaires (distorsion harmonique, fréquence dépendant de l'amplitude) qui seront étudiés en détail dans ce cours. Le Chapitre 2 détaille les techniques qui serviront à l'analyse des systèmes vibratoires non linéaires effectuées dans les chapitres qui suivent. Le forçage harmonique de l'oscillateur de Duffing est étudié dans le Chapitre 3. La résonance primaire et les résonances secondaires (super-harmonique et sous-harmonique) sont étudiées par la méthode des échelles multiples. Le chapitre 4 est consacré à l'analyse de systèmes auto-oscillants. Le cas d'un amortissement négatif induit par le frottement est étudié par la méthode des échelles multiples et celui de l'instabilité aéroélastique d'une aile par la méthode de moyennisation de Bogoliubov-Krylov.

**Notions abordées.** Oscillateur de Duffing (résonances primaire et secondaires), stabilité linéaire des systèmes dynamiques, oscillations auto-entretenues (induites par le frottement, flottement), méthodes perturbatives (échelles de temps multiples, moyennage de Bogoliubov-Krylov).

<footnote n="1">1. FISE ≡ Formation Initiale sous Statut Étudiant.</footnote>

<footnote n="2">2. GSI ≡ Génie des Systèmes Industriels.</footnote>

<footnote n="3">3. IMC ≡ Ingénieries Mécanique et Conception.</footnote>
## Vibrations des structures (2017-aujourd'hui)

J'ai repris et fait évoluer ce cours depuis 2017. Il était initialement donné par Jean-Mathieu Mencik (Professeur de Université à l'INSA Centre Val de Loire). Il comporte 12h de Cours Magistraux, 9h20 de Travaux Dirigés et 8h de Travaux Pratiques pour les étudiants en 4e année du département GSI. Je n'interviens pas directement en Travaux Pratique.

L'objectif du cours est de former les étudiants aux méthodes d'analyses et de calcul en dynamique des structures. En particulier, après un premier chapitre rappelant la notion d'oscillateur harmonique, l'analyse modale des systèmes discrets à plusieurs degrés de liberté et des systèmes continus simples est abordée. La réponse de ces systèmes à des sollicitations harmoniques est également étudiée.

**Notions abordées.** Oscillations amorties des systèmes à plusieurs degrés de liberté, Oscillations des systèmes continus (barre en traction-compression et poutre en flexion), forçage harmonique, réponse en fréquence, analyse modale.

## Matériaux pour l'ingénieur (2015-aujourd'hui)

Ce cours est également une création. Il comporte 9h20 de Cours Magistral et 20h de Travaux Dirigés pour les étudiants en 4e année du département GSI. Je prenais initialement en charge l'intégralité du cours mais depuis deux ans je n'interviens plus que la première partie consacrée à la physique des matériaux métalliques. Les reste est dispensé par mes collègues Arnaud Frachon, Michaël Caliez et Michel Grattons (Maîtres de Conférences à l'INSA Centre Val de Loire).

L'objectif du cours est d'étudier un large spectre de la physique des matériaux métalliques mais également le comportement des matériaux des polymères. Le cours est divisée en quatre partie : (i) physique des matériaux métalliques, (ii) comportement des matériaux métalliques, (iii) comportement des polymères et (iv) mécanique du collage.

**Notions abordées.** Physique des matériaux métalliques : composition et structure des matériaux, diffusion à l'état solide, thermochimie, propriétés mécaniques. Élasticité et plasticité 3D des métaux, fatigue des métaux, nature chimique et recyclage des polymères, mise en œuvre, caractérisation et comportement des polymères, dimensionnement d'un assemblage collé.

## Qualité et maîtrise statistique des processus (2015-2022).

J'ai repris et fait évoluer ce cours depuis 2015. Initialement donné par Marie-Laure Gobert (Maîtresse de Conférences à l'INSA Centre Val de Loire), il comporte 9h de Cours Magistraux et 8h20 de Travaux Dirigés pour les étudiants en 3e année FISE et 4e année FISA 4 du département GSI.

L'objectif du cours est (i) de sensibiliser les étudiants à la notion de qualité en milieu industriel dans le cadre des normes ISO 9000 et ISO 9001, (ii) d'introduire la notion de capabilité et (iii) de présenter un outil d'amélioration de la qualité particulièrement adapté aux processus de production de moyennes et grandes séries : la maîtrise statistique des processus (MSP).

**Notions abordées.** Capabilité, cartes de contrôle aux moyennes (carte de Shewhart, carte de type Shewhart, carte EWMA), carte de contrôle à la dispersion, contrôle aux attributs.

<footnote n="4">4. FISE ≡ Formation Initiales sous Statut Apprenti.</footnote>
## Mécanique du point (2017-aujourd'hui)

Ce cours comporte 10h40 de Cour Magistral, 24h de Travaux Dirigés et 5h20 de Travaux Pratiques pour les étudiants en 1re année du département STPI<fnref n="5" />. Les Cours Magistraux sont dispensés par Frédéric Mabilat (Professeur Agrégé à l'INSA Centre Val de Loire) et Marie-Laure Gobert. J'interviens dans les Travaux Dirigés et les Travaux Pratiques.

Il s'agit d'un cours classique de mécanique du point pour des étudiants-ingénieurs débutant leur cursus universitaire. Le cours est divisé en trois parties :

1. Dynamique et énergétique en coordonnées cartésiennes : Principe Fondamental de la dynamique (PFD), Théorème de l'Énergie Cinétique (TEC).

2. Dynamique (PFD) et énergétique (TEC) en coordonnées cylindriques ; Théorème du moment cinétique ; statique des solides.

3. Energie potentielle et mécanique (théorème de l'énergie mécanique) ; Force de rappel élastique exercée par un ressort ; Oscillations mécaniques

**Notions abordées.** Initiation au calcul vectoriel, dynamique du point, cinématique et dynamique en coordonnées cylindriques, forces conservatives et énergie, énergie mécanique et théorèmes associés, oscillations créées par un ressort.

## Conception assistée par ordinateur (2015-aujourd'hui)

Ce cours comporte 8h40 de Travaux Dirigés et 15h20 de Travaux Pratiques pour les étudiants en 1re année du département STPI. Il est sous la responsabilité de Flavien Baudou (Professeur Agrégé à l'INSA Centre Val de Loire). J'interviens exclusivement en Travaux Pratiques.

L'objectif du cours est l'apprentissage d'un logiciel de C.A.O. 3D (Solidworks).

**Notions abordées.** Conception (par extrusion, révolution ou balayage) des différentes pièces d'un assemblage simple, réalisation de l'assemblage, réalisation d'une mise en plan côtée, utilisation d'une bibliothèque de composants.

## 6.3 À L'UNIVERSITÉ DE TOURS

### Fluid-structure interactions (2015-2018)

Ce cours fut une création pour les étudiants en 2e année du Master International « Electronics and Mechanical Engineering » et était dispensé en anglais. Il comportait 6h de Cours Magistraux et 6h de Travaux Dirigés.

Ce cours considère le problème vibroacoustique du rayonnement acoustique de structures élastiques couplées à un fluide. L'objectif est de mettre en évidence dans les cas de structures simples comme les plaques, les lois essentielles qui régissent les transferts vibroacoustiques. Le problème est abordé de façon analytique.

**Notions abordées.** Équations de l'acoustique linéaire, équations d'onde dans les plaques, couplage vibra-acoustique (cas d'une plaque infinie) dans les cas de fluides légers et de fluides lourds, courbes de dispersion.

<footnote n="5">5. STPI ≡ Sciences et Technologies Pour l'Ingénieur.</footnote>
## 6.4 À L'UNIVERSITÉ DU MANS

Au cours de mes trois années de doctorat à l'Université de Mans et dans le cadre de mon contrat doctoral, j'ai eu l'opportunité d'effectuer une mission complémentaire « enseignement ». Ce fut pour moi l'occasion d'enseigner à plusieurs niveaux dans les formations proposées à l'UFR Sciences et Techniques de l'Université du Mans. Ces enseignements sont brièvement décrits ci-dessous.

### Onde plane (2013-2015)

Responsable du module : Jean-Pierre Dalmont (Professeur des Universités à l'Université de Mans).

Volume horaire : 26h de Cours et de Travaux dirigés pour des étudiants et 2ᵉ année de Licence STI<fnref n="6" />.

**Contenu scientifique.** Les équations de l'acoustique linéaire (en onde plane) sont établies à partir des lois fondamentales de la physique : principe fondamentale de la dynamique (i.e. équation d'Euler), la conservation de la masse et une équation thermodynamique (i.e. équation d'état). Les solutions des équations de l'acoustique sont étudiées dans le domaine temporel puis fréquentiel. Sont ensuite étudiés successivement les notions de réflexion et transmission des ondes, de mode propre et de résonance, d'énergie, d'intensité et de puissance des ondes acoustiques. Enfin les mécanismes de dissipation et leur prise en compte dans les équations de l'acoustique sont abordés.

### Notions élémentaires d'acoustique (2013-2015)

Responsable du module : Jean-Pierre Dalmont.

Volume horaire : 24h de Travaux Pratiques pour des étudiants et 1ᵉ année de Licence STI.

**Contenu scientifique.** Corde vibrante, tube de Kundt, réflexion et transmission des ondes, interférences, propagation d'ondes ultrasonores, pendules couplés, sources à directivité variable (monopôle, dipôle, quadripôle), propagation guidée avec écoulement.

### Perception du son et codage (2013-2015)

Responsable du module : André Almeida (Maître de Conférences à l'Université de Mans).

Volume horaire : 12h de Travaux Pratiques pour des étudiants et 1ᵉ année de Licence STI.

**Contenu scientifique.** Notion de son musical, sources cohérentes/incohérentes et battement, initiation à l'audiologie, échelle de sonie.

### Électroacoustique avancée (2013)

Responsable du module : Bruno Gazengel (Professeur des Universités à l'Université de Mans).

Volume horaire : 24h de Travaux Pratiques pour des étudiants et 2ᵉ année du DEUST Vibration, Acoustique et Signal (VAS).

**Contenu scientifique.** Mesure de l'impédance acoustique d'un guide d'onde, étude des caractéristiques d'un silencieux acoustique, étude des charges d'un haut-parleur, filtrage électrique d'un haut-parleur, directivité de système électroacoustique, caractérisation de microphones de mesure.

<footnote n="6">6. SPI ≡ Sciences Pour l'Ingénieur.</footnote>
# Physique des instruments de musique (2013-2015)

Responsable du module : Jean-Pierre Dalmont.

Volume horaire : 12h de Travaux Pratiques pour des étudiants et 1re année du Master Ingénierie Mécanique et Acoustique (IMA).

**Contenu scientifique.** Mesure et analyse des courbes d'impédance d'entrée de différents instruments à vent, oscillations auto-entretenues (exemple d'un modèle simplifié de clarinette), couplage structure-cavité (exemple du couplage table-caisse d'une guitare classique).

**Contribution personnelle notable.** Élaboration du TP « Oscillations auto-entretenues : modèle simplifié de clarinette ».
BIBLIOGRAPHIE PERSONNELLE

ARTICLES DANS DES REVUES INTERNATIONALES À COMITÉ DE LECTURE

[a1] A. ALMEIDA, B. BERGEOT, C. VERGEZ et B. GAZENGEL. « Analytical Determination of the Attack Transient in a Clarinet With Time-Varying Blowing Pressure ». *Acta Acustica united with Acustica* 101.5 (2015), p. 1026-1038. DOI: 10.3813/AAA.918897.

[a2] B. BERGEOT. « Scaling law for the slow flow of an unstable mechanical system coupled to a nonlinear energy sink ». *Journal of Sound and Vibration* 503 (2021), p. 116109. ISSN: 0022-460X. DOI: doi.org/10.1016/j.jsv.2021.116109.

[a3] B. BERGEOT. « Effect of stochastic forcing on the dynamic behavior of a self-sustained oscillator coupled to a non-linear energy sink ». *International Journal of Non-Linear Mechanics* 150 (2023), p. 104351. DOI: 10.1016/j.ijnonlinmec.2023.104351.

[a4] B. BERGEOT, A. ALMEIDA, B. GAZENGEL, C. VERGEZ et D. FERRAND. « Response of an artificially blown clarinet to different blowing pressure profiles ». *The Journal of the Acoustical Society of America* 135.1 (2014), p. 479. DOI: 10.1121/1.4835755.

[a5] B. BERGEOT, A. ALMEIDA, C. VERGEZ et B. GAZENGEL. « Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure ». *Nonlinear Dynamics* 73.1-2 (2013), p. 521-534. DOI: 10.1007/s11071-013-0806-y.

[a6] B. BERGEOT, A. ALMEIDA, C. VERGEZ et B. GAZENGEL. « Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure : influence of noise ». *Nonlinear Dynamics* 74.3 (2013), p. 591-605. DOI: 10.1007/s11071-013-0991-8.

[a7] B. BERGEOT et S. BELLIZZI. « Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network ». *Nonlinear Dynamics* 94.2 (2018), p. 1501-1522. DOI: 10.1007/s11071-018-4438-0.

[a8] B. BERGEOT et S. BELLIZZI. « Steady-state regimes prediction of a multi-degree-of-freedom unstable dynamical system coupled to a set of nonlinear energy sinks ». *Mechanical Systems and Signal Processing* 131 (2019), p. 728-750. DOI: https://doi.org/10.1016/j.ymssp.2019.05.045.

[a9] B. BERGEOT, S. BELLIZZI et S. BERGER. « Dynamic behavior analysis of a mechanical system with two unstable modes coupled to a single nonlinear energy sink ». *Communications in Nonlinear Science and Numerical Simulation* 95 (2021), p. 105623. ISSN: 1007-5704. DOI: https://doi.org/10.1016/j.cnsns.2020.105623.

[a10] B. BERGEOT, S. BELLIZZI et B. COCHELIN. « Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment ». *International Journal of Non-Linear Mechanics* 78 (2016), p. 72-89. DOI: 10.1016/j.ijnonlinmec.2015.10.006.

[a11] B. BERGEOT, S. BELLIZZI et B. COCHELIN. « Passive suppression of helicopter ground resonance instability by means of a strongly nonlinear absorber ». *Advances in Aircraft and Spacecraft Science* 3.3 (2016), p. 271-298.

[a12] B. BERGEOT, S. BELLIZZI et B. COCHELIN. « Passive suppression of helicopter ground resonance using nonlinear energy sinks attached on the helicopter blades ». *Journal of Sound and Vibration* 392 (2017), p. 41-55. DOI: 10.1016/j.jsv.2016.12.039.

[a13] B. BERGEOT et S. BERGER. « Fast-slow analysis of passive mitigation of self-sustained oscillations by means of a bistable nonlinear energy sink ». *Physica D : Nonlinear Phenomena* 460 (2024), p. 134063. DOI: 10.1016/j.physd.2024.134063. URL: https://hal.science/hal-04280444.

[a14] B. BERGEOT, S. BERGER et S. BELLIZZI. « Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis ». *Journal of Vibration and Control* 24.15 (2017), p. 3487-3511. ISSN: 1077-5463. DOI: 10.1177/1077546317707101.