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

Christophe Vergez

Analytical prediction of delayed Hopf bifurcations in a simplified stochastic model of reed musical instruments

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

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.

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 th 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 dependance)

A(y)

Intermediate quantity used throughout the paper and defined by Eq. [START_REF] Stratonovich | Topics In the Theory of Random Noise[END_REF]. 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 ŷdyn stoch,a )

ν, σ Magnitude of the stochastic forcing (σ = ν/ √ 2)

ϵ Rate of linear increasing of the bifurcation parameter α i , ω i , F i Modal parameters of the i th mode (damping, eigenfrequency, modal factor respectively)

ζ
Embouchure parameter (physical model of clarinet)

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 [START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF]. In this "static" context, periodic solutions and their stability can be also determined using for example the harmonic balance method [START_REF] Fréour | Numerical continuation of a physical model of brass instruments: Application to trumpet comparisons[END_REF][START_REF] Colinot | Multiple two-step oscillation regimes produced by the alto saxophone[END_REF] or orthogonal collocation [START_REF] Terrien | Flute-like musical instruments: A toy model investigated through numerical continuation[END_REF][START_REF] Gilbert | Minimal blowing pressure allowing periodic oscillations in a simplified reed musical instrument model: Bouasse-Benade prescription assessed through numerical continuation[END_REF]. 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 timevarying 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 [START_REF] Bergeot | Response of an artificially blown clarinet to different blowing pressure profiles[END_REF]: 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 clarinet1 [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF][START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: Influence of noise[END_REF]. 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 [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF][START_REF] Baesens | Noise effect on dynamic bifurcations: The case of a period-doubling cascade[END_REF][START_REF] Jansons | Stochastic calculus: Application to dynamic bifurcations and threshold crossings[END_REF][START_REF] Berglund | Noise-induced phenomena in slow-fast dynamical systems[END_REF]. 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 impor-tant 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 [START_REF] Silva | MoReeSC: A framework for the simulation and analysis of sound production in reed and brass instruments[END_REF]. 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 [START_REF] Chaigne | Acoustics of Musical Instruments[END_REF]. 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? 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 2. Section 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 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 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 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 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 6.

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

Physical problem

Dimensionless equations, Eq. ( 11)

Amplitude-Phase representation, Eq. ( 18)

Stochastic slow dynamic, Eq. ( 29)

Analysis of the dynamic bifurcation, Sect. 2 Single reed instrument model with a white noise forcing

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 [START_REF] Benade | Fundamentals of musical acoustics[END_REF][START_REF] Fletcher | The physics of musical instruments[END_REF][START_REF] Chaigne | Acoustics of Musical Instruments[END_REF]: the cane reed and the air-column inside the instrument. While blowing air through the reed channel into the instrument, the musician provides a quasistatic source of energy. The instrument and the player constitute an autonomous dynamical system. When the triv-ial equilibrium solution of this system becomes unstable, a sound is produced [START_REF] Wilson | Operating modes of the clarinet[END_REF][START_REF] Fletcher | Nonlinear theory of musical wind instruments[END_REF][START_REF] Silva | Interaction of reed and acoustic resonator in clarinet-like systems[END_REF].

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 [START_REF] Backus | Small-vibration theory of the clarinet[END_REF][START_REF] Ollivier | Idealized models of reed woodwinds. part 1 : Analogy with bowed string[END_REF]. 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 (ω) = Z(ω)U (ω), ( 1
)
where ω 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(ω). For computational reasons, the series is truncated to N modes, where N is an integer:

Z(ω) = N n=1 F n jω ω 2 n + jωα n ω n -ω 2 , (2)
with F n , ω n and α 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. ( 2) can be written in the time domain:

pn + α n ω n ṗn + ω 2 n p n = F n u, ∀n ∈ [1, N ], (3)
with u the inverse Fourier transform of U and p n is such that p = N i=1 p n , where p is the inverse Fourier transform of P [START_REF] Silva | Interaction of reed and acoustic resonator in clarinet-like systems[END_REF] 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 [START_REF] Hirschberg | Aero-acoustics of wind instruments[END_REF][START_REF] Dalmont | Nonlinear characteristics of single-reed instruments: Quasistatic volume flow and reed opening measurements[END_REF]. A polynomial expansion of this relation is often written in the neighborhood of the equilibrium solution (i.e. the mean flow) [START_REF] Kergomard | Elementary considerations on reedinstrument oscillations[END_REF]:

u(t) = u eq + c 1 p(t) + c 2 p(t) 2 + c 3 p(t) 3 , ( 4
)
with , where γ is the dimensionless pressure in the mouth of the musician and ζ a dimensionless parameter accounting for many embouchure parameters. By Eq. ( 4), Eq. ( 3) can be written using only the pressure p as follows

u eq = ζ(1 -γ) √ γ the mean volume flow, c 1 = ζ 3γ-1 2γ 1 2 , c 2 = -ζ 3γ+1
pn + α n ω n ṗn + ω 2 n p n + F n ṗf (p, γ) = 0, ∀n ∈ [1, N ], (5)
where f (p, γ) = -∂u ∂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 (5) becomes

p + α 1 ω 1 ṗ + ω 2 1 p + F 1 ṗf (p, γ) = 0 (6)
In this case (N = 1), since p 1 = p, p 1 is replaced by p in Eq. ( 6). 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 γ. The stability of the trivial equilibrium solution p = 0 with respect to γ is classically analyzed by looking at the sign of the eigenvalues real parts of the Jacobian matrix of Eq. ( 6) written in the state-space form. This leads to the following expression for the static Hopf bifurcation point

γst = 1 3 + 2α 1 ω 1 α 1 ω 1 + α 2 1 ω 2 1 + 3ζ 2 F 2 1 9ζ 2 F 2 1 (7)
corresponding to the value of γ for which the two complex conjugate eigenvalues become with positive real parts. Note that in the lossless case (i.e. α 1 = 0) the static bifurcation parameter is γst = 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. (7) are given for example by Kergomard et al. [START_REF] Kergomard | Calculation of the spectrum of self-sustained oscillators using a variable truncation method: Application to cylindrical reed instruments[END_REF] and Silva et al. [START_REF] Silva | Interaction of reed and acoustic resonator in clarinetlike systems[END_REF].

Compared to the model given by [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: Influence of noise[END_REF], this article studies the case of a linearly increasing blowing pressure γ, with additional stochastic excitation as detailed in Sect. 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 found for example in [START_REF] Øksendal | Stochastic Differential Equations: An Introduction with Applications[END_REF].

The one-mode model described by Eq. ( 6) is now considered with a linearly increasing blowing pressure γ = εt+γ 0 . The variation of the blowing pressure is assumed to be small during a period T 1 = 2π ω1 which is ensured by 0 < ε ≪ 1. Because one assumes also that ω 1 ≫ 1 (without loss a generality in musical context) the time variation of the mouth pressure γ is neglected in the time derivative of the mean flow [START_REF] Benade | Fundamentals of musical acoustics[END_REF]. Moreover, the model is now forced by a stochastic excitation νξ t where ν is the noise level and ξ t is assumed to be a unitary idealized white noise process with zero mean, i.e.

E [ξ t ] = 0 and E [ξ t ξ t+τ ] = δ(τ ) ( 8
)
where δ is the Dirac delta function and E [{.}] denotes the ensemble average2 of {.}. Finally, the one-mode stochastic model is written as follows

pt + α 1 ω 1 ṗt + ω 2 1 p t + F 1 ṗt f (p t , γ t ) = νξ t (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. ( 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. ( 9). First, a new bifurcation parameter y t = γ t -γst is considered to obtain a system whose trivial solution becomes unstable at y = 0 (i.e. ŷst = 0). Then the time rescaling t → t ′ = ω 1 t is introduced. Recalling that a normalized white noise ξ t is defined as the time derivative of the normalized Wiener process W t and using the scaling property of Wiener process (see [START_REF] Øksendal | Stochastic Differential Equations: An Introduction with Applications[END_REF], Chap. 2) which states that W t ω 1 and 1

√ ω1 W t are equivalent, we have

ξ t = dW t dt ∼ ω 1 dW t ′ ω 1 dt ′ ∼ √ ω 1 dW t ′ dt ′ = √ ω 1 ξ t ′ . ( 10
)
Therefore, by using also {} for the derivation with respect to t ′ and denoting t ′ by t for the sake of conciseness, Eq. ( 9) takes the form of the following stochastically forced selfexcited oscillator

pt + h (p t , ṗt , y t ) + p t = νξ t ( 11
)
where

h (p t , ṗt , y t ) = α 1 ṗt + F 1 ω 1 ṗt f p t , y t + γst (12)
and ν =

ν ω 3/2 1
. The time evolution of y t is given by

y t = ϵt + y 0 ( 13
)
where y 0 = γ 0 -γst and ϵ = ε ω1 with 0 < ϵ ≪ 1 (because ω 1 ≫ 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 ω 1 t.

Following Roberts and Spanos [START_REF] Roberts | Stochastic averaging: An approximate method of solving random vibration problems[END_REF] the stochastic averaging method [START_REF] Stratonovich | Topics In the Theory of Random Noise[END_REF][START_REF] Khas'minskii | A Limit Theorem for the Solutions of Differential Equations with Random Right-Hand Sides[END_REF], whose general formulation is recalled in Appendix B, is used to obtain the stochastic slow dynamics of Eq. [START_REF] Colinot | Numerical simulation of woodwind dynamics: investigating nonlinear sound production behavior in saxophone-like instruments[END_REF]. To apply the method, the response process (p t , ṗ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 + φ t ) (14a) ṗt = -x t sin (t + φ t ) (14b)
The desired form of Eq. (14b) requires that

ẋt cos ϕ t -x t φt sin ϕ t = 0, (15)
with

ϕ t = t + φ t which yields φt = ẋt x t cos ϕ t sin ϕ t and ẋt = x t φt sin ϕ t cos ϕ t . ( 16
)
Then differentiation of Eq. (14b) leads to

pt = -x t cos ϕ t -ẋt sin ϕ t -x t φ cos ϕ t . (17)
Finally, the substitution of ( 14) and ( 17) into [START_REF] Colinot | Numerical simulation of woodwind dynamics: investigating nonlinear sound production behavior in saxophone-like instruments[END_REF] and the use of ( 16) yields ẋt = h (x t cos ϕ t , -x t sin ϕ t , y t ) sin ϕ t -νξ t sin ϕ t (18a)

φt = h (x t cos ϕ t , -x t sin ϕ t , y t ) cos ϕ t x t -νξ t cos ϕ t x t (18b)
which has a similar form as (70) in Appendix B with

x t = (x t , φ t ) T , η t = (ξ t , ξ t ) T , f (x t , t) =   h (x t cos ϕ t , -x t sin ϕ t , y t ) sin ϕ t 1 x t h (x t cos ϕ t , -x t sin ϕ t , y t ) cos ϕ t   , g(x t , t) =   -ν sin ϕ t 0 0 -ν cos ϕ t x t   (19)
with () T the transpose operator. Following the approach described in Appendix B, the drift vector m (see Eq. ( 72)) is first computed. The first term T av {f } = (F (x t , y t ), G(x t , y t ))

T , with

F (x t , y t ) = F 1 ω 1 ζ 3(y t + γst ) -1 4(y t + γst ) 1/2 - α 1 2 x t + F 1 ω 1 ζ 3 (γ st + y t + 1) 128 (γ st + y t ) 5/2 x 3 t (20) G(x t , y t ) = 0, (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 ϕ t are slow variables with respect to the unitary eigenfrequency of the considered dimensionless simplified reed instrument model [START_REF] Colinot | Numerical simulation of woodwind dynamics: investigating nonlinear sound production behavior in saxophone-like instruments[END_REF] which therefore corresponds to a period equal to 2π. Then, the second part of the drift vector m (see again Eq. ( 72)) is determined as

T av 0 -∞ E ∂(gη) ∂x t (gη) t+τ dτ = 1 2π 2π 0 0 -∞ E ∂(gη) ∂x t (gη) t+τ dτ dt = 1 2π 0 -∞ 2π 0 δ(τ ) ν 2 cos ϕt cos ϕt+τ xt -ν 2 sin ϕt cos ϕt+τ +cos ϕt sin ϕt+τ x 2 t dtdτ = ν 2 4xt 0 ( 22
)
in which the symbol E [] disappeared due to definitions [START_REF] Berglund | Pathwise description of dynamic pitchfork bifurcations with additive noise[END_REF] and one used the fact that

0 -∞ δ(τ ) cos(τ )dτ = 1 2 +∞ -∞ δ(τ ) cos(τ )dτ = 1 2 .
The final expression of the drift vector is therefore

m =   F (x t , y t ) + ν 2 4x t 0   . ( 23
)
Now the expression of the diffusion matrix σ is determined from (73) starting by computing

T av +∞ -∞ E (gη) t (gη) T t+τ dτ = 1 2π 2π 0 +∞ -∞ E (gη) t (gη) T t+τ dτ dt = +∞ -∞ ν 2 δ(τ ) cos(τ ) 2 -ν 2 δ(τ ) sin(τ ) 2xt ν 2 δ(τ ) sin(τ ) 2xt ν 2 δ(τ ) cos(τ ) 2x 2 t dτ =    ν 2 2 0 0 ν 2 2x 2 t    . ( 24
)
in which again definitions [START_REF] Berglund | Pathwise description of dynamic pitchfork bifurcations with additive noise[END_REF] have been used. From (24) a possible solution of (73) is

σ =    ν √ 2 0 0 1 x t ν √ 2    . ( 25
)
Therefore, according to (71) the equations governing the stochastic slow dynamics of [START_REF] Colinot | Numerical simulation of woodwind dynamics: investigating nonlinear sound production behavior in saxophone-like instruments[END_REF] are

dx t = F (x t , y t ) + ν 2 4x t dt + ν √ 2 dW t (26a) dφ t = 1 x t ν √ 2 dW t . ( 26b
)
in which the time evolution of the amplitude x t is uncoupled from that of the phase φ t .

Assuming a weak noise level (i.e. ν ≪ 1) and therefore neglecting the term ν 2 4xt in (26a) and denoting σ = ν √ 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 + σdW t . ( 27
)
The scaling and time shift properties of Wiener process (see [START_REF] Øksendal | Stochastic Differential Equations: An Introduction with Applications[END_REF], Chaps. 7) allow us to write that

1 √ ϵ (W ϵt+y0 -W y0 ) ∼ 1 √ ϵ W ϵt ∼ W t and therefore 1 √ ϵ dW y ∼ dW t . (28)
Finally, from ( 27) and ( 28) we obtain the final form of the slow dynamics as

dx y = 1 ϵ F (x y , y)dy + σ √ ϵ dW y . ( 29
)
where x t and y t are now denoted x y and y respectively to emphasize that Eq. ( 29) is a non-autonomous onedimensional Itô stochastic differential equation with respect to x y depending on y which acts as the time variable.

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 (11) becomes a supercritical pitchfork bifurcation for Eq. ( 27). Therefore, if the averaged system ( 27) is a good approximation of ( 11), the dynamic Hopf bifurcation point of Eq. ( 11) is very close to the dynamic pitchfork bifurcation point of [START_REF] Ollivier | Idealized models of reed woodwinds. part 1 : Analogy with bowed string[END_REF]. The analytical expression of the latter is obtained in this section in both deterministic and stochastic cases.

As stated by Stocks et al. [START_REF] Stocks | Influence of random fluctuations on delayed bifurcations: The case of additive white noise[END_REF], the dynamic bifurcation point 3 has two main definitions: it is (1) the value of the bifurcation parameter y when E x 2 y crosses a predefined threshold x 2 th (here the initial value is chosen) or (2) the value of the bifurcation parameter y at which, on average, the random variable x 2 y crosses the same given threshold x 2 th . 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. [START_REF] Baesens | Noise effect on dynamic bifurcations: The case of a period-doubling cascade[END_REF] and [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: Influence of noise[END_REF] for a clarinet-like system). Definition 4.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 ŷdyn , is defined as the value of y t for which E x 2 y 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.

Solution of the linearized averaged system

The linearized version of (29) with respect to x y around 0 is considered

dx y = 1 ϵ a(y)x y dy + σ √ ϵ dW y ( 30
)
where, from Eq. ( 20), one has

a(y) = ∂F ∂x y (0, y) = F 1 ω 1 ζ 3(y + γst ) -1 4(y + γst ) 1/2 - α 1 2 . ( 31
)
In order to solve [START_REF] Silva | Interaction of reed and acoustic resonator in clarinetlike systems[END_REF], the deterministic differential equation associated to the latter is first considered as

dx y dy = 1 ϵ a(y)x y ( 32
)
whose solution, denoted x det (y), is

x det (y) = x y0 e 1 ϵ (A(y)-A(y0)) (33)
where

A(y) = ζF 1 (γ st + y -1) γst + y 2ω 1 - α 1 2 y. ( 34
)
is the antiderivative of a(y). Following a well-known method for solving stochastic differential equations, we apply the Itô formula (see Eq. (69) in Appendix A) to the function f (x y , y) = x y e -1 ϵ (A(y)-A(y0)) which is constant for the deterministic Eq. (32) (indeed through Eq. ( 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 ∂f ∂y + 1 ϵ a(y)x y ∂f ∂x vanishes, we obtain

df (x y , y) = σ √ ϵ e -1 ϵ (A(y)-A(y0)) dW y ( 35
)
and integrating [START_REF] Taillard | An Analytical Prediction of the Bifurcation Scheme of a Clarinet-Like Instrument: Effects of Resonator Losses[END_REF] from y 0 to y yields the following solution of Eq. ( 30)

x y = x det (y) + σ √ ϵ e 1 ϵ A(y) y y0 e -1 ϵ A(y ′ ) dW y ′ . ( 36
)
Note that the first and second terms of the right-hand side of ( 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 E gdW y = 0 assuming some properties for the function g (see [START_REF] Øksendal | Stochastic Differential Equations: An Introduction with Applications[END_REF], Chap. 3) that are respected by e -1 ϵ A(y) . Therefore we have

E y y0 e -1 ϵ A(y) dW y = 0. ( 37
)
That means that the expected value of x y is the solution of the associated deterministic equation [START_REF] Spiegel | Mathematical Handbook of Formulas and Tables[END_REF], i.e. E [x y ] = x det (y). Now E x 2 y is computed from [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF], that yields

E x 2 y = E (x det (y)) 2 + E σ 2 ϵ e 2 ϵ A(y) y y0 e -1 ϵ A(y ′ ) dW y ′ 2 + E 2 x det (y) σ √ ϵ e 1 ϵ A(y) y y0 e -1 ϵ A(y ′ ) dW y ′ . ( 38
)
Using again the property of the Itô integral which states that E gdW y = 0, the third term in the right-hand side of ( 38) is equal to zero. The second term is simplified using the Itô isometry (see Corollary 3.1.7 in [START_REF] Øksendal | Stochastic Differential Equations: An Introduction with Applications[END_REF], Chap. 3) which states that E gdW y 2 = E g 2 dy assuming again some properties for the function g that are resected by e -1 ϵ A(y) . Therefore we obtain the final expression of the expected value of the squared amplitude

E x 2 y = D(y) + S(y) (39)
where

D(y) = (x det (y)) 2 (40a) S(y) = σ 2 ϵ e 2 ϵ A(y) y y0 e -2 ϵ A(y) dy (40b)
are respectively the deterministic and the stochastic parts of the expected value of the squared amplitude. From (39) and according to [START_REF] Berglund | Pathwise description of dynamic pitchfork bifurcations with additive noise[END_REF], three regimes can be distinguished for a system such as (29) undergoing a dynamic pitchfork bifurcation:

• Regime I: S(y) ≪ 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 (32) which undergoes the greatest possible bifurcation delay.

• Regime II: S(y) ≫ D(y) with a noise level not too high. The noise can no longer be neglected but it remains small enough for linearized model [START_REF] Silva | Interaction of reed and acoustic resonator in clarinetlike systems[END_REF] to remain valid. A bifurcation delay still exists but it is reduced compared to Regime I.

• Regime III: S(y) ≫ 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 [START_REF] Berglund | Pathwise description of dynamic pitchfork bifurcations with additive noise[END_REF], the notion of bifurcation delay becomes meaningless.

The domains of existence of each regime are given explicitly in Sect. 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.

Dynamic bifurcation point in the deterministic case

In this section we assume that the system is in Regime I and therefore E x 2 y = D(y) = (x det (y)) 2 (with x det (y) given by Eq. ( 33)). Consequently, from Definition 4.1, the deterministic dynamic bifurcation point, denoted ŷdyn det , is a solution of x y0 = x y0 e 1 ϵ (A(y)-A(y0)) and therefore of

A(y) = A(y 0 ). ( 41
)
From [START_REF] Stratonovich | Topics In the Theory of Random Noise[END_REF], Eq. ( 41) is solved (details are given in Appendix C). It is shown that, in addition to the trivial solution y = y 0 , Eq. (41) has two other solutions. The one that corresponds to the deterministic dynamic bifurcation point is

ŷdyn det = X 2 2 -γst , ( 42
)
where X 2 is given by Eq. (78). We will see in Fig. 2 that the maximum of ŷdyn det is obtained for y 0 = -γ st and in this case Eq. (42) reduces to

ŷdyn det = α 1 ω 1 + α 2 1 ω 2 1 + 4ζ 2 F 2 1 2 4ζ 2 F 2 1 -γst (43)
It is important to be aware that even if ŷdyn det has a finite value for y 0 = -γ st , the function F (x, y) in Eq. ( 29) diverge at y = -γ st . Therefore, Eq. (43) must be understood as a limit value for y 0 very close to -γ st (from above). However, for values of y very close to -γ st , the linearization of Eq. ( 29) with respect to x (leading to Eq. ( 30)) is valid if x is very close to zero. In other words, Eq. ( 43) is valid to predict the deterministic bifurcation point of a numerical simulation of Eq. ( 29) only if y 0 and x y0 are very close 42), as a function of the initial condition y0. The opposite of the static bifurcation point γst , given by Eq. ( 7), is depicted by a vertical dashed line. The red points correspond to the initial conditions used in Fig. 3. The set of parameters (45) is used.

to -γ st and 0 respectively. This corresponds to the situation where the musician begins to blow gently into the instrument.

In the lossless case, i.e. α 1 = 0 and γst = 1 3 , Eq. ( 42) becomes

ŷdyn det = 1 2 1 -y 0 -1 + y 0 (2 -3y 0 ) (44)
which depends only on the initial value y 0 . Note first that Eq. (44) yields ŷdyn det = 2/3 if y 0 = -γ st = -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 γ = γst + ŷdyn det = 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. α 1 > 0).

Then, in cases with and without damping, ŷdyn det does not depend on the slope ϵ which may seem counterintuitive.

The deterministic dynamic bifurcation point ŷdyn det , given by Eq. (42), is plotted in Fig. 2 as a function of the initial condition y 0 with the following set of parameters:

ϵ = 0.002, ω 1 = 1000 rad • s -1 , α 1 = 0.02, F 1 = 1200 s -1 and ζ = 0.2. ( 45
)
The figure shows that ŷdyn det is a decreasing function with respect to y 0 and, as previously mentioned, the maximum is obtained for y 0 = -γ st . Moreover, the red points correspond to the initial conditions used in Fig. 3. In the latter, numerical simulations of the deterministic differential deterministic differential equation associated to Eq. ( 29) are shown with a logarithmic scale for the vertical axis. The same parameters are as in Fig. 2 and the initial values are: x(y 0 ) = 0.01 and y 0 = -0.36, -0.32, . . . , -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 ŷ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 ŷdyn det . The green and red points are used respectively to highlight the considered initial conditions and the corresponding dynamic bifurcation points predicted by Eq. ( 42). This shows that the theoretical results presented in Fig. 2 predict the deterministic bifurcation points measured on numerical simulations.

In the lossless case the expression of the deterministic dy-

namic bifurcation point ŷdyn

det is given by (44) which depends only on the initial value y 0 . The deterministic bifurcation point ŷdyn det given by Eq. ( 42) is plotted as a function of the initial condition y 0 for five values of the damping coefficient in Fig. 4 and as a function of the parameter ζ for four values of the damping coefficient α 1 and for three values of the initial condition y 0 in Fig. 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 ζ and for initial conditions y 0 close to the opposite of the static bifurcation point γst (see Fig. 5). In general one has 0.1 < ζ < 0.4 for a clarinet, 0.25 < ζ < 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. 42) as a function of the parameter ζ for four of the damping coefficient (α1 = 0, 0.01, 0.02 and 0.05), for three value of the initial condition (y0 = -γ st , -3γ st /4 and -γ st /2) with ω1 = 1000 rad•s -1 and F1 = 1200 s -1 .

Dynamic bifurcation point in the stochastic case

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

x 2 y0 = σ 2 ϵ e 2 ϵ A(y) y y0 e -2 ϵ A(y) dy. ( 46
)
The first step is to obtain the approximate expression of the integral

I(y 0 , y) = y y0 e -2 ϵ A(y) dy. ( 47
)
For this purpose A(y) (given by Eq. ( 34)) is expanded in a second order Taylor series around 0 (the static bifurcation point), that leads to

A(y) ≈ A(0) + y dA dy (0) + y 2 2 d 2 A dy 2 (0) (48a) = A(0) + ya(0) + y 2 2 a ′ (0) (48b)
where the second term in the right-hand side of (48b) vanishes by definition4 and a ′ (0) = da dy (0). The approximate expression of I(y 0 , y) is therefore given at order 2 by

I(y 0 , y) = y y0 e -2 ϵ A(0)+ y 2 2 a ′ (0) dy = 1 2 πϵ a ′ (0) e -2 ϵ A(0) × erf y a ′ (0) ϵ -erf y 0 a ′ (0) ϵ ( 49
)
where erf is the error function.

Since we consider 0 < ϵ ≪ 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 a ′ (0) ϵ ≪ -1 ⇒ erf y 0 a ′ (0) ϵ ≈ -1 (50)
and for y > 0 (we are interested in the dynamic bifurcation point which is by definition larger than the static bifurcation point)

y a ′ (0) ϵ ≫ 1 ⇒ erf y a ′ (0) ϵ ≈ 1. (51)
Therefore, from Eqs. ( 50) and (51), Eq. ( 49) reduces to

I(y 0 , y) = πϵ a ′ (0) e -2 ϵ A(0) . ( 52
)
Eq. (52) highlights that I(y 0 , y) is now independent of y and y 0 and simply denoted I. Therefore Eq. (46) becomes

A(y) = K, ( 53
) with K = -ϵ ln σ - ϵ 2 ln I ϵ + ϵ ln x y0 (54)
and the function A(y) given by [START_REF] Stratonovich | Topics In the Theory of Random Noise[END_REF].

Eq. ( 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 ( 55
)
with

a 1 = ζ 2 F 2 1 4ω 2 1 , a 2 = 1 4 (3γ st -2)ζ 2 F 2 1 ω 2 1 -α 2 1 , a 3 = (γ st -1)(3γ st -1)ζ 2 F 2 1 4ω 2 1 -α 1 K, a 4 = (γ st -1) 2 γst ζ 2 F 2 1 4ω 2 1 -K 2 (56)
and γst given by Eq. ( 7) (details are given in Appendix D). Eq. ( 55) is solved analytically using Cardano's formula (see again Appendix D for details). One obtains ŷdyn stoch,a = r 0 if σ < σ 3 , (57a)

r 2 if σ 3 < σ < σ 2 (57b)
where r 0 and r 2 are given by Eq. ( 85), as depicted in Fig. 13. The expressions of σ 2 and σ 3 are given respectively by Eqs. (82) to (84).

A second approximate expression of the stochastic dynamic bifurcation point is obtained using (48b) in Eq. ( 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

ŷdyn stoch,b = 2 - (γ st ) 3/2 ω 1 ϵ ζF 1 (3γ st + 1) × 3 ln(γ st ) 2 + ln 8πω 1 ϵζF 1 (3γ st + 1) + 4 ln(σ) -4 ln (x y0 ) 1/2 . ( 58
)
Note that both in Eqs. ( 57) and ( 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 ϵ.

Deterministic and stochastic dynamic bifurcation points, given by Eqs. ( 42), ( 57) and (58) respectively, are plotted in Fig. 6 as functions of the noise level σ. The deterministic bifurcation point ŷdyn det is plotted for two values of the initial condition y 0 , i.e. y 0 = -γ st and -3γ st /4. The stochastic bifurcation points are plotted for two values of the parameter ϵ, i.e. ϵ = 0.002 and 0.01. The other parameters are given by (45) with in addition x y0 = 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 (48b) is valid. Secondly, Fig. 6 shows graphically the domains of existence of the regimes listed at the end of the Sect. 4.1.

The following definitions are chosen for the boundary values (with respect to the noise level σ) between the different regimes previously mentioned: the boundary value between Regime I and Regime II, denoted σ I/II , corresponds to the intersection between ŷdyn det and ŷdyn stoch,a . This means in Fig. 6 that if σ < σ I/II (Regime I), ŷdyn det is the actual dynamic bifurcation point. On the contrary, if σ > σ I/II (Regime II), ŷdyn stoch,a or ŷdyn stoch,b (depending on the approximation retained) is the actual dynamic bifurcation point.

Domains of existence of the regimes

The boundary value between Regime I and Regime II, with respect to the noise level σ and denoted σ I/II , is obtained solving

D(y) = S(y). ( 59
)
with respect to σ. Therefore, from (40) and ( 52), we find

σ I/II = x y0 a ′ (0)ϵ π 1/4 e 2 ϵ (A(0)-A(y0)) . ( 60
)
The second boundary value between Regime II and Regime III, again with respect to the noise level σ denoted σ II/III , corresponds to the value of σ for which ŷdyn stoch,a and ŷdyn stoch,b do not exist anymore5 (see Fig. 6). The expression of σ II/III is obtained by noting that ŷdyn stoch,a and 61) respectively, corresponding to situations plotted in Fig. 6. with γst given by Eq. ( 7). Using Eqs. ( 60) and ( 61) respectively the values of σ I/II and σ II/III corresponding to situations plotted in Fig. 6 are given in Tab. 1.

y 0 = -γ st y 0 = -3γ st /
By means of Eqs. ( 60) and ( 61) one can also depict the regions of existence of each of the three regimes in the (y 0 , ϵ, σ)-space (see Fig. 7(a)). A representation of these regions is also given in the (y 0 , σ)-plane for ϵ = 0.002 and in the plane (ϵ, σ)-plane for y 0 = -0.34 in Figs. 7(b) and 7(c) respectively. One can see that the boundary between Regime II and Regime III mainly depends on σ. On the contrary, the boundary between Regime I and Regime II depends simultaneously on σ, y 0 and ϵ. That means that the use of Eq. ( 42) (the deterministic bifurcation point) or Eq. ( 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. 7 shows that the dependence on ϵ is the most important. Indeed, the above value σ I/II = 3.15•10 -10 has been obtained for ϵ 0.002. The value becomes σ I/II = 7.21 • 10 -17 for ϵ = 10 -3 and σ I/II = 1.44 • 10 -69 for ϵ = 2 • 10 -4 . In other words, if ϵ = 2 • 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 [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: Influence of noise[END_REF].

Probability density function of the stochastic averaged amplitude

In addition to being able to obtain the expected value of the squared amplitude (39), Eq. ( 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. ( 29). Indeed, a known result of stochastic calculus (see e.g. [START_REF] Klebaner | Introduction to stochastic calculus with applications[END_REF], Chap. 4) is that the Itô integral of a deterministic function f (t), i.e. t 0 f (t)dW t , is a Gaussian process with mean equal to zero and variance v = t 0 f (t) 2 dt. Therefore, from (36), x y is a Gaussian process with mean E [x y ] = x det (y) and variance

E x 2 y -E [x y ] 2 = S(y) = σ 2 ϵ e 2 ϵ A(y) y y0 e -2 ϵ A(y) dy (62)
with y y0 e -2 ϵ A(y) dy given by (52). The associated Gaussian PDF is therefore

ρ(x, y) = 1 2πS (y) e -(x-x det (y)) 2 2S(y) . (63)
This means that, for a given value of y, the probability that x y lies between x and x + dx is ρ(x, y)dx. Then, we define the function R(y) = ρ(x y0 , y) for which the probability that x y lies between x y0 and x y0 + dx is R(y)dx. From (63), the expression of the function R(y) is

R(y) = 1 2πS(y) e -(xy 0 -x det (y)) 2 2S(y) . ( 64
)
In the Regime I we have D(y) = x det (y) 2 ≫ S(y) and the function R(y) can be approximated by a Dirac delta function as R(y) ≈ δ(x y0 -x det (y)). Denoting ℓ(y) = x y0x det (y) the roots of which are y 0 and ŷdyn det (see Eq. ( 41)) and using the composition property of Dirac delta function one obtains

R(y) ≈ R I (y) = δ(y -y 0 ) |ℓ ′ (y 0 )| + δ(y -ŷdyn det ) |ℓ ′ (ŷ dyn det )| (65)
when the system evolves according to Regime I. In the Regime II, because D(y) ≪ S(y), Eq. (64) becomes

R(y) ≈ R II (y) = 1 2πS(y) e - x 2 y 0 2S(y) . ( 66
)
We have therefore which vanishes for x y0 = S(y) (it can be shown that this corresponds to a maximum of R ′ II (y)) and therefore, according to the Definition 4.1, for y = ŷdyn stoch,a . In Fig. 8 the function R(y) is plotted using parameters (45) and three values of the noise level, i.e. σ = 10 -4 , σ = 10 -6 and σ = 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. 10 in the next section). In Figs. 8(a) and 8(b) one can see that the maximum of R(y) is obtained for y = ŷdyn stoch,a as predicted by Eq. (67). Eq. ( 65) is illustrated in Fig. 8(c) in which we can see that in Regime I the function R(y) becomes a Dirac delta function translated to y = ŷdyn det . 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. (65) does not appear.

R ′ II (y) = S ′ (y) x 2 y0 -S(y) √ 2πS(y) 5/2 e -

Comparison between theoretical results and numerical simulation

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 ŷdyn det of the deterministic dynamic bifurcation point given by Eq. ( 42) and those of the stochastic bifurcation points, ŷdyn stoch,a and ŷdyn stoch,b , given by Eqs. ( 57) and ( 58) respectively.

To achieve that, after getting the time series of p t from the numerical integration of (11) 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 divided into N intervals and on each of them the maximum of p t (max(p) i with i = 1, . . . , N ) and the mean of y t (mean(y) i with i = 1, . . . , 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 E A 2 y is computed over the 50 realizations. Simultaneously, the expected value E x 2 y of x y is computed over 50 realizations of the numerical integration of the averaged equation [START_REF] Silva | Interaction of reed and acoustic resonator in clarinet-like systems[END_REF].

The results are shown in Fig. 9 which represents: (on the left) a plot of E A 2 y (dashed red line), E x 2 y (blue line), the static bifurcation diagram (red) computed in Appendix E and ŷdyn det , ŷdyn stoch,a and ŷdyn stoch,b 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 E x 2 y = x y0 and using a logarithm scale for the E x 2 y -axis. The set of parameters (45) is used and, from top to bottom: σ = 10 -4 , σ = 10 -6 and σ = 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. 9(a) and 9(b) the system is in Regime II. Therefore, the noise cannot be ignored and the dynamic bifurcation must be determined by ŷdyn stoch,a or ŷdyn stoch,b . The figure shows that the latter provides a good approximation of the dynamic bifurcation point. Indeed, through the Definition 4.1, the dynamic bifurcation point is the abscissa of the point of intersection between E [x y ] and the horizontal at E [x y ] = x y0 . The logarithmic scale used in Fig. 9 (right column) allows us to locate this intersection and to see that the vertical green line, corresponding to ŷdyn stoch,a , passes also through this intersection and the vertical magenta dashed line, corresponding to ŷdyn stoch,b , passes a little on the left (all the more as σ is decreased). Of course the expression ŷdyn stoch,a is more accurate than ŷdyn stoch,b because it corresponds to a lesser degree of approximation. In Fig. 9(c) the system is in Regime I, consequently, the noise can be neglected and the dynamic bifurcation point is determined by ŷdyn det . 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 E [x y ] and the horizontal line at E [x y ] = x y0 .

The position in the (y 0 , σ)-plane of the three situations depicted in Fig. 9 is shown in Fig. 10 and superimposed to the regions of existence of the Regimes I, II and III.

Comparison in term of probability density function

In this section the theoretical PDF ρ(x, y) (see Eq.( 63)) is compared to two histograms: one in built from the full order system [START_REF] Colinot | Numerical simulation of woodwind dynamics: investigating nonlinear sound production behavior in saxophone-like instruments[END_REF] and the other from the stochastic slow dynamics [START_REF] Silva | Interaction of reed and acoustic resonator in clarinet-like systems[END_REF]. 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 ỹ (chosen arbitrary to be ỹ = 0.15). The comparison, presented in Fig. 11, is performed for three values of the noise level (the same as chosen previously in Figs. 8 and9) and shows en excellent agreement between theoretical and numerical results. In Fig. 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 ỹ 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.

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 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 [START_REF] Colinot | Numerical simulation of woodwind dynamics: investigating nonlinear sound production behavior in saxophone-like instruments[END_REF]). 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 [START_REF] Chaigne | Acoustics of Musical Instruments[END_REF]. (68)

where m and σ are real functions and W t is the so-called Wiener process.

Let f (x t , t) ∈ C 2 (R, R + ) (i.e. f is twice continuously differentiable on (R, R + )). Then f (x t , t) is also Itô process, and

df (x t , t) = ∂f ∂t dt + ∂f ∂x dx t + 1 2 ∂ 2 f ∂x 2 σ(x t , t) 2 dt (69a) = ∂f ∂t + m(x t , t) ∂f ∂x + σ(x t , t) 2 2 ∂ 2 f ∂x 2 dt + σ(x t , t) ∂f ∂x dW t . ( 69b
)
Eq. ( 69) is the 1-dimensional Itô's formula (for more details and proof see [START_REF] Øksendal | Stochastic Differential Equations: An Introduction with Applications[END_REF], Chap. 4).

B General formulation of the stochastic averaging method

In this appendix the stochastic averaging method [START_REF] Stratonovich | Topics In the Theory of Random Noise[END_REF][START_REF] Khas'minskii | A Limit Theorem for the Solutions of Differential Equations with Random Right-Hand Sides[END_REF] is briefly described. For that we consider the following system of differential equations in standard form ẋt = f (x t , t) + g(x t , t)η t (70) where x t ∈ R n . If the deterministic vector function f (x t , t) ∈ R n and matrix function g(x t , t) ∈ R n ×R n satisfy certain requirements [START_REF] Khas'minskii | A Limit Theorem for the Solutions of Differential Equations with Random Right-Hand Sides[END_REF] and the elements of the vector η t are broadband processes, with zero means, then the slow (or averaged) dynamics of Eq. ( 70) may be approximated by the following Itô equations

dx t = m(x t )dt + σ(x t , t)dW t , (71)
where W t ∈ R n is a vector of n Wiener processes. The vector m and the matrix σ are called drift vector and diffusion matrix respectively and defined by

m = T av f + 0 -∞ E ∂(gη) ∂x t (gη) t+τ dτ , ( 72
)
where ∂(gη) ∂x is the Jacobian matrix of gη, and

σσ T = T av +∞ -∞ E (gη) t (gη) T t+τ dτ ( 73
)
where {.} T and E [{.}] denotes respectively the transpose and the expected value of {.}. T av is an averaging operator defined as follows

T av {.} = lim T →+∞ 1 T t0+T t0 {.} dt. ( 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 av becomes a classical Krylov-Bogolyubov time averaging over one period T 0 , i.e.

T av {.} = 1 T 0 t0+T0 t0 {.} dt (75)
and the result is independent of t 0 .

C Derivation of the expression of the deterministic dynamic bifurcation point

In this appendix we give the details of the résolution of Eq. ( 41). First we state y + γst = X 2 and then it can be shown that, using Eq. ( 34), Eq. ( 41) takes the following form

(X -X 0 ) ζF 1 X 2 + (ζF 1 X 0 -α 1 ω 1 )X -ζF 1 + ζF 1 X 2 0 -α 1 X 0 ω 1 = 0 (76)
where X 2 0 = y 0 + γst . Obviously X = X 0 and therefore y = y 0 is a solution of (76). The second term of the product in the left-hand side of (76) is a second order polynomial equations whose roots X 1 and X 2 are

X 1 = 1 2ζF 1 α 1 ω 1 -ζF 1 X 0 -α 2 1 ω 2 1 + 2α 1 ζF 1 X 0 ω 1 + ζ 2 F 2 1 (4 -3X 2 0 ) (77)
and

X 2 = 1 2ζF 1 α 1 ω 1 -ζF 1 X 0 + α 2 1 ω 2 1 + 2α 1 ζF 1 X 0 ω 1 + ζ 2 F 2 1 (4 -3X 2 0 ) . ( 78
)
The initial value y 0 is always chosen to be larger than -γ st (because the mouth pressure γ 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

ŷdyn det = X 2 2 -γst . ( 79
)

D Derivation of the expression of the stochastic dynamic bifurcation point

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

K = ζF 1 γst + y 2ω 1 y -1 + γst - α 1 y 2 . ( 80
)
To obtain a solvable cubic form (i.e. without square root), Eq. ( 80) is transformed into

K + α 1 y 2 2 = ζ 2 F 2 1 (γ st + y) 4ω 2 1 y -1 + γst 2 (81)
which yields Eq. ( 55).

In the remaining of this appendix the Cardano's method (see e.g. [START_REF] Spiegel | Mathematical Handbook of Formulas and Tables[END_REF]) is used to solve the latter, i.e. a 1 y 3 + a 2 y 2 + a 3 y + a 4 = 0.

First, the following parameters are introduced

p = - a 2 2 3a 2 1 + c 3 a 1 and q = a 2 27a 1 2a 2 2 a 2 1 - 9a 3 a 1 + a 4 a 1 .
The discriminant ∆ is defined as ∆ = -4p 3 + 27q 2 . Then:

1. If ∆ < 0, one root is real and two are complex conjugate.

2. If ∆ = 0, all roots are real and at least two are equal.

3. If ∆ > 0, all roots are real and unequal.

A typical example of the discriminant ∆, plotted as a function of the noise level σ, is depicted in Fig. 12 for a typical set of parameters. It can be shown that ∆ can be expressed as a fourth order polynomial equation with respect to ln σ which has two distinct roots ln (σ 1 ) = 2α with k = 0, 1, 2. (85) Fig. 13 shows the roots r k (k = 0, 1, 2) as functions of the noise level σ using again the parameters (45) and x y0 = 0.01. The stochastic dynamic bifurcation point, denoted ŷdyn stoch,a , is equal to r 0 if σ < σ 3 and to r 2 if σ 3 < σ < σ 2 . This choice is justified by means of a numerical resolution which shows that this is the unique positive solution of Eq. ( 80).

3 1 ω 3 1 -27α 1 γst ζ 2 F 2 1 ω 1 54ζ 2 F 2 1 ω 1 ϵ - 2 α 2 1 ω 2 1 + 3ζ 2 F 2 1 3/2 -9α 1 ζ 2 F 2 1 ω 1 54ζ 2 F 2 1 ω 1 ϵ - 1 2

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. ( 6) 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. ( 27) is considered without noise et with a constant bifurcation parameter y, i.e.

dx t dt = F (x t , y), ( 86
)
where the function F (x t , y t ) is given by Eq. ( 20). The fixed points x e of (86) are obtained by solving F (x, y) = 0. We obtain three solutions: the trivial solution x e 1 = 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 e 2 , its expression is (87)

where the expression of ŷst is given by Eq. ( 7). In the lossless case with α 1 = 0 and γst = 1 3 , Eq. (87) reduces to x e 2 = 4 √ 2(3y + 1) y 9y+12 . The trivial solution corresponds to zero equilibrium position of (6) whereas the non trivial solution characterizes its periodic steady-state regimes.

As we know, the trivial fixed point x e 1 is stable if y < 0 and unstable if y > 0. One can be shown that the non trivial fixed point x e 2 exists only for y > 0 and is stable.

γ, γ t

  Dimensionless blowing pressure (the subscript () t indicates that γ varies with time) γst Value of γ corresponding to a static Hopf bifurcation y t Time-varying dimensionless blowing overpressure with respect to the static Hopf bifurcation (= γ t -γst ) ŷdyn det 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 ŷdyn stoch,a Analytical prediction of the value of y at which the system undergoes a dynamic Hopf bifurca-1 tion in the stochastic case (the subscript () a means first level of approximation) ŷdyn stoch,b

  8γ

3 2 and c 3 = -ζ γ+1 16γ 5 2

 3232

Figure 2 .

 2 Figure 2. Deterministic dynamic bifurcation point ŷdyndet , given by Eq. (42), as a function of the initial condition y0. The opposite of the static bifurcation point γst , given by Eq. (7), is depicted by a vertical dashed line. The red points correspond to the initial conditions used in Fig.3. The set of parameters (45) is used.

1 Figure 3 .

 13 Figure 3. Numerical simulations of the deterministic differential equation associated to Eq. (29) using a logarithm scale for the vertical axis. The set of parameters (45) is used with, in addition, x(y0) = 0.01 (the latter is depicted by a horizontal dashed line) and y0 = -0.36, -0.32, . . . , -0.16. The static bifurcation point ŷst = 0 is depicted by a vertical gray dashed line. The dynamic bifurcation point ŷdyn det is defined for a given initial condition by the value of y for which |x(y)| = |x(y0)|, 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. (42).

Figure 4 .

 4 Figure 4. The deterministic bifurcation point ŷdyn det given by Eq. (42) as a function of the initial condition y0 for five values of the damping coefficient, i.e. α1 = 0, 0.025 . . . , 0.1 with ω1 = 1000 rad•s -1 , F1 = 1200 s -1 and ζ = 0.2. Red dashed lines indicate the value of -γ st for each value of αn.

Figure 5 .

 5 Figure 5. The deterministic bifurcation point ŷdyndet given by Eq. (42) as a function of the parameter ζ for four of the damping coefficient (α1 = 0, 0.01, 0.02 and 0.05), for three value of the initial condition (y0 = -γ st , -3γ st /4 and -γ st /2) with ω1 = 1000 rad•s -1 and F1 = 1200 s -1 .

Figure 6 .

 6 Figure 6. Deterministic and stochastic dynamic bifurcation points, given by Eqs. (42), (57) and (58) respectively, as functions of the noise level σ. The deterministic bifurcation point ŷdyn det is plotted for two values of the initial condition y0, i.e. y0 = -γ st and -3γ st /4. The stochastic bifurcation points are plotted for two values of the parameter ϵ, i.e. ϵ = 0.002 and 0.01. The other parameters are given by (45) with in addition xy 0 = 0.01.

Figure 7 .

 7 Figure 7. Regions of existence of the Regimes I, II and III (a) in the (y0, ϵ, σ)-space, (b) in the (y0, σ)-plane for ϵ = 0.002 and (c) in the plane (ϵ, σ)-plane for y0 = -0.34. The expressions of σ I/II and σ II/III are given by Eqs. (60) and (61) respectively.

Figure 8 .

 8 Figure 8. The function R(y) defined by Eq. (64). The quantities ŷdyn det , ŷdyn stoch,a and ŷdyn stoch,b are depicted by vertical green line and vertical red dashes line respectively. The parameters (45) are used and (a) σ = 10 -4 , (b) σ = 10 -6 and (c) σ = 10 -11 . Moreover y0 = -0.34, xy 0 = 0.01 and ϵ = 0.002.

Figure 9 . 1 Figure 10 .

 9110 Figure 9. On the left: a plot with the expected values of A 2 y (red dashed line) and x 2 y (blue line) obtained using 50 realizations of the numerical integration of Eqs (11) and (29) respectively, the stable branches of the static bifurcation diagram (SBD, gray line) computed in Appendix E and ŷdyn det , ŷdyn stoch,a and ŷdyn stoch,b 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 E [x 2 y ] = xy 0 and using a logarithm scale for the E [x 2 y ]-axis. The set of parameters (45) is used and, from (a) to (c), σ = 10 -4 , σ = 10 -6 and σ = 10 -11 . Moreover y0 = -0.34, xy 0 = 0.01 and ϵ = 0.002.

Figure 11 .

 11 Figure 11. Comparison between theoretical PDF ρ(x, y) given by Eq. (63) and two histograms: one in built from the full order system[START_REF] Colinot | Numerical simulation of woodwind dynamics: investigating nonlinear sound production behavior in saxophone-like instruments[END_REF] and the other from the stochastic slow dynamics[START_REF] Silva | Interaction of reed and acoustic resonator in clarinet-like systems[END_REF]. In both cases, 2000 realizations of the system are computed and we take, for each realization, the value of the considered random variable (xy or pt) for a given value of y denoted ỹ (chosen arbitrary to be ỹ = 0.15). The comparison is performed for three values of the noise level: (a) σ = 10 -4 , (b) σ = 10 -6 and (c) σ = 10 -11 . Moreover y0 = -0.34, xy 0 = 0.01 and ϵ = 0.002.

ln   2 √ 2 √Figure 12 .Figure 13 .

 221213 Figure 12. Discriminant ∆ of (55) as a function the natural logarithm of the noise level ln σ. The set of parameters (45) is used.

  st +y+1) ω1(γ st +y)5/2

Table 1 .

 1 Values

σ I/II and σ II/III , computed through Eqs. (60) and (

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[START_REF] Dalmont | Oscillation and extinction thresholds of the clarinet: Comparison of analytical results and experiments[END_REF][START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF][START_REF] Taillard | An Analytical Prediction of the Bifurcation Scheme of a Clarinet-Like Instrument: Effects of Resonator Losses[END_REF].

We recall that the ensemble average consists in repeating the same measurement many times, and in calculating the average over them.

Equations governing the stochastic slow dynamics

Sometimes called exit value in the literature.

Indeed, at the bifurcation one has a(0) = ∂F ∂x (0, 0) = 0.

This corresponds to the fact that the solutions of Eq. (53) for both cubic and quadratic approximations become complex for σ > σ II/III .

All numerical simulations of the Itô stochastic differential equations are performed using the function ItoProcess of the Wolfram Mathematica software.

Code availability On request.

A The 1-dimensional Itô's formula

Let the following Itô differential equation

Declarations

Funding Not applicable.

Conflict of interest

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

Availability of data and material On request.