S Terrien

C Vergez

B Fabre

D Barton

Calculation of the steady-state oscillations of a flute model using the orthogonal collocation method

In this paper we exploit the method of numerical continuation combined with orthogonal collocation to determine the steady-state oscillations of a model of flute-like instruments that is formulated as a nonlinear neutral delay differential equation. The delay term in the model causes additional complications in the analysis of the behaviour of the model, in contrast to models of other wind instruments which are formulated as ordinary differential equations. Fortunately, numerical continuation provides bifurcation diagrams that show branches of stable and unstable static and periodic solutions of the model and their connections at bifurcations, thus enabling an in depth analysis of the global dynamics. Furthermore, it allows us to predict the thresholds of the different registers of the flute-like instrument and thus to explain the classical phenomenon of register change and the associated hysteresis.

Introduction

Flute-like musical instruments present complex behaviour, some aspects of which remain poorly understood. For instance, one can cite the fact that a regime change can lead either to a periodic or to a quasiperiodic solution, the generation mechanism of such quasiperiodic regimes, or the way in which some parameters affect the regime change thresholds (and thus the importance of the associated hysteresis). The study of a physical model of this instrument class provides the opportunity to explore the influence of different parameters and, by comparison with experimental results, to improve the knowledge of the physical mechanisms of sound production.

Many studies have dealt with the calculation of self-sustained oscillations of flute models. Timedomain simulations (see for example [START_REF] J W Coltman | Time-domain simulation of the flute[END_REF][START_REF] J W Coltman | Jet offset, harmonic content, and warble in the flute[END_REF][START_REF] Auvray | Regime change and oscillation thresholds in recorderlike instruments[END_REF]) allow the computation of both transients and steady-state periodic and non-periodic oscillations, and are thus informative about the diversity of dynamical regimes. Nevertheless, these methods are very sensitive to initial conditions, and do not necessary provide information about coexistence of multiple solutions. Some analytical methods were proposed by Schumacher and Fletcher [START_REF] Schumacher | Self-sustained oscillations of organ flue pipes: an integral equation solution[END_REF][START_REF] Fletcher | Sound production by organ flue pipes[END_REF], but require drastic simplifications of the model. Moreover, due to both their complexity and the limitation to a few number of harmonics, these methods are not suited for a systematic study of the instrument behaviour.

Numerical methods for the calculation of periodic solutions of nonlinear dynamical systems have been developed for decades in the field of nonlinear differential equations: they include in particular the harmonic balance method (HBM), based on a frequency-domain discretisation of the solution (see for example [START_REF] Nakhla | A piecewise harmonic balance technique for determination of periodic response of nonlinear systems. Circuits and Systems[END_REF]), and the orthogonal collocation method [START_REF] Doedel | Lecture notes on numerical analysis of nonlinear equations[END_REF][START_REF] Doedel | AUTO: A program for automatic bifurcation analysis of autonomous systems[END_REF], based on a time-domain discretisation.

Since they can be combined with numerical continuation methods [START_REF] Krauskopf | Numerical continuation methods for dynamical systems[END_REF], these techniques are particularly interesting. Indeed, the HBM or the orthogonal collocation considered alone compute periodic solutions x for different values of the model parameters λ (see figure 1-(a)), each point being independent from the others. Numerical continuation methods rely on the implicit function theorem (see for example [START_REF] Krauskopf | Numerical continuation methods for dynamical systems[END_REF]) to say that these points lay on a continuous solution branch. Therefore, they can follow the evolution of a given The computation of the different branches leads to bifurcation diagrams, which represent static and periodic solutions of the model, as functions of a parameter of interest (the so-called continuation parameter). Since bifurcation diagrams give information about the coexistence of multiple solutions and the presence of unstable solutions, such diagrams provide a more global knowlegde of the system dynamics.

For the first time in 1989, the harmonic balance method was introduced in musical acoustics by Gilbert et al. [START_REF] Gilbert | Calculation of the steady-state oscillations of a clarinet using the harmonic balance technique[END_REF]. This technique has then been used in various works focusing on single reed musical instruments [START_REF] Fritz | Some aspects of the harmonic balance method applied to the clarinet[END_REF][START_REF] Farner | Contribution to harmonic balance calculations of self-sustained periodic oscillations with focus on single-reed instruments[END_REF][START_REF] Dalmont | Reed instruments, from small to large amplitude periodic oscillations and the helmholtz motion analogy[END_REF] and brass instruments [START_REF] Menguy | Weakly nonlinear gas oscillations in air-filled tubes; solutions and experiments[END_REF][START_REF] Gilbert | A simulation tool for brassiness studies[END_REF]. Recent studies [START_REF] Karkar | Toward the systematic investigation of periodic solutions in single reed woodwind instruments[END_REF][START_REF] Karkar | Oscillation threshold of a clarinet model: a numerical continuation approach[END_REF] have then demonstrated the benefits of combining this technique with numerical continuation methods to understand the behaviour of clarinet-like instruments. In parallel, software packages have been developed (among which [START_REF] Doedel | AUTO: A program for automatic bifurcation analysis of autonomous systems[END_REF][START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF][START_REF] Engelborghs | DDE Biftool: a Matlab package for bifurcation analysis of delay differential equations[END_REF][START_REF] Dhooge | Matcont: a matlab package for numerical bifurcation analysis of odes[END_REF]), that combine computation of periodic solutions and numerical continuation methods.

Compared to models of other wind instruments, the state-of-the-art model for flute-like instruments contains an additional delay term due to the presence of an air jet in the excitation mechanism, and is moreover a neutral delayed dynamical system (that is to say it includes a delay in the highest derivative -see section 2). Such systems are mathematically more complex than ordinary differential equations governing, for example, the behaviour of reed instruments. In principle, the harmonic balance method can be used to solve this kind of systems [START_REF] Junca | Interaction between periodic elastic waves and two contact nonlinearities[END_REF]. However, the current work does not aim to develop from scratch a new software package specifically dedicated to the flute model, but rather to take advantage of existing validated numerical tools. Nevertheless, the different numerical tools cited above do not allow calculation and continuation of periodic solutions of neutral equations, and the associated numerical schemes must be adapted.

A first attempt to study flute-like instruments through numerical continuation, using a toymodel (mathematically a more simple non-neutral delayed dynamical system), has highlighted that such an approach provides new kind of information, explaining behaviours which can be related to some experimental phenomena [START_REF] Terrien | Flutelike musical instruments: A toy model investigated through numerical continuation[END_REF]. It thus suggests to study the state-of-the-art model of flute-like instruments in the same way, and therefore to introduce numerical tools adapted to neutral systems. This is done thanks to the contribution of Barton [START_REF] Barton | Collocation schemes for periodic solutions of neutral delay differential equations[END_REF], whose work allows the analysis of neutral systems within the DDE-Biftool environment.

We first recall in section 2 the equations of the state-of-the-art physical model of flute-like instruments. We then describe in section 3 the orthogonal collocation method and the predictorcorrector approach, which respectively allow the computation and the continuation of periodic solutions of neutral delay dynamical systems. Finally, an application to the model of flute-like instruments highlights in section 4 the valuable contribution of such an approach to the understanding of the phenomenon of "register change" and its associated hysteresis, two typical features of flute-like instruments.

2 State-of-the-art model for flutes

Modeling of the mechanism of sound production

As for other wind instruments, the establishment of auto-oscillation in flute-like instruments results from the coupling between an exciter and a resonator constituted by the air column contained in the pipe. The peculiarity of flute-like instruments lies in the nature of the exciter, whose a schematic representation is provided in figure 2.

Whereas the excitation mechanism of other wind instruments involves the vibration of a solid element (such as a cane reed or the musician's lips), it is constituted here by the nonlinear interaction of an air jet with a sharp edge called "labium". More precisely, when the musician blows in the instrument, a jet is created at the channel exit. As highlighted in figure 2, this channel is a part of the instrument in the case of recorders or flue organ pipes, and formed by the player's lips for transverse flutes. As the jet is naturally unstable, any perturbation is amplified while convected along the jet, from the channel exit to the labium. The jet-labium interaction then makes the jet oscillate around the labium, leading to an alternative flow injection inside and outside the instrument, that constitutes an aeroacoustic pressure source for the resonator. The acoustic waves created in the resonator disrupt back the air jet at the channel exit, thus sustaining the auto-oscillation of the instrument. A simplified representation of the jet behaviour is provided in figure 3. As reviewed by Fabre et al. [START_REF] Fabre | Physical modeling of flue instruments: A review of lumped models[END_REF], this mechanism of sound production can be represented by a feedback loop system, represented in figure 4, which involves a separate description of the main physical phenomena. We briefly recall here the physical mechanism related to each element and the associated equations.

Receptivity of the jet

Once the auto-oscillations established, the perturbation of the jet is provided by the acoustic waves in the pipe. The reaction of the jet in terms of flow transverse displacement η 0 at the channel exit is called "receptivity" (see figure 3):

η 0 (t) = h U j v ac (t), ( 1
)
where h is the height of the channel (highlighted in figure 2), U j the jet velocity, and v ac (t) the acoustic velocity at the resonator inlet.

Figure 4: Basic modeling of sound production mechanism in flute-like instruments, as a system with a feedback loop [START_REF] Chaigne | Acoustique des instruments de musique (Acoustics of musical instruments)[END_REF][START_REF] Fabre | Physical modeling of flue instruments: A review of lumped models[END_REF].

Amplification and convection of the perturbations along the jet

As the jet is naturally unstable, the initial perturbation η 0 is amplified during its convection from the channel exit to the labium. Experimental study by de la Cuadra [START_REF] De | The sound of oscillating air jets: Physics, modeling and simulation in flute-like instruments[END_REF] has shown that the exponential amplification of the perturbation with the convection distance W , initially described by Rayleigh [START_REF] Rayleigh | The theory of sound second edition[END_REF] in the case of small perturbations of an infinite jet, remains a reasonable approximation for flute-like instruments:

η(W, t) = η 0 (t -τ )e α i W , ( 2
)
where W is here the distance between the channel exit and the labium (see figure 2), and α i an amplification coefficient, for which de la Cuadra [START_REF] De | The sound of oscillating air jets: Physics, modeling and simulation in flute-like instruments[END_REF] determined the following empirical expression: α i ≈ 0.4 h . The delay τ , introduced by the duration of convection of the initial perturbation η 0 along the jet, is related both to the distance W and to the convection velocity c v of transversal perturbations on the jet: τ = W cv . Both theoretical and experimental results [START_REF] Rayleigh | The theory of sound second edition[END_REF][START_REF] De La Cuadra | Visualization and analysis of jet oscillation under transverse acoustic perturbation[END_REF][START_REF] Nolle | Sinuous instability of a planar jet: propagation parameters and acoustic excitation[END_REF] have shown that c v is related to the jet velocity itself through: c v ≈ 0.4U j .

Aeroacoustic source: jet-labium interaction

The jet-drive model, initially proposed by Coltman [START_REF] J W Coltman | Jet drive mechanisms in edge tones and organ pipes[END_REF] and followed by Verge [START_REF] Verge | Jet oscillations and jet drive in recorder-like instruments[END_REF][START_REF] Verge | Sound production in recorder-like instruments. ii. a simulation model[END_REF], is based on a representation of the source term as a force. More precisely, the jet oscillates around the labium, which splits the related air flow in a part Q in which enters in the pipe, and a part Q out which comes out of the pipe, as highlighted in figure 3. These two flow sources in phase opposition are considered to be localised at points close to the labium and separated by a distance δ d evaluated by Verge in [START_REF] Verge | Jet oscillations and jet drive in recorder-like instruments[END_REF] as:

δ d ≈ 4 π √ 2hW .
As this distance remains small compared to the acoustic wavelength, the air between the two flow injection points is assumed to be incompressible [START_REF] Chaigne | Acoustique des instruments de musique (Acoustics of musical instruments)[END_REF]. This mass of air acts as a dipolar pressure source ∆p src (t) on the air column contained in the pipe. The second Newton's law leads to:

∆p src (t) = - ρδ d W H • dQ in dt , (3)
where H is the width of the excitation window, and ρ the air density. The flow Q in , related to the transversal deflection of the jet η(W, t) at the labium, can be calculated from the velocity profile of the air jet U (y):

Q in = H y 0 -η(W,t) -∞ U (y)dy, (4)
where y 0 is the offset between the labium position and the jet centerline (see figure 2). Following Segouffin [START_REF] Ségoufin | Experimental study of the influence of the mouth geometry on sound production in a recorder-like instrument: Windway length and chamfers[END_REF], we assume that a Bickley profile of half width b = 2h 5 is a correct approximation of the jet velocity profile U (y) at the labium:

U (y) = U j sech 2 y b . (5)
Assuming that the centerline jet velocity U j varies in a quasi-static way, one finally obtains the following modeling of the pressure source:

∆p src (t) = ρδ d bU j W • d dt tanh η(W, t) -y 0 b . (6)

Vortex shedding at the labium : nonlinear losses

Between the channel exit and the labium of flue instruments, the presence of an important transversal flow induced by the acoustic field in the pipe can cause a flow separation at the labium, resulting in vortex shedding [START_REF] Fabre | Vortex shedding in steady oscillation of a flue organ pipe[END_REF]. This phenomenon corresponds to an energy dissipation, and can be represented by an additional (nonlinear) pressure variation ∆p los near the labium:

∆p los (t) = - ρ 2 v ac (t) α vc 2 sgn(v ac (t)) (7)
where α vc is a vena contracta factor, estimated at 0.6 in the case of a sharp edge, and sgn represents the sign function.

Finally, the pressure source can be written as:

∆p(t) = ∆p src (t) + ∆p los (t). (8)

Passive response of the resonator

The acoustical response of the air column contained in the pipe, excited by the pressure source described above, is represented through the input admittance Y in (ω) = V ac (ω)/∆P (ω), where ω is the pulsation, and V ac (ω) and ∆P (ω) are respectively the acoustic velocity at the pipe inlet and the pressure source, written in the frequency domain. The input admittance Y in (ω) is represented as a sum of resonance modes, including a mode at zero frequency (the so-called uniform mode [START_REF] Chaigne | Acoustique des instruments de musique (Acoustics of musical instruments)[END_REF]):

Y in (ω) = a 0 b 0 jω + c 0 + p-1 k=1 a k jω ω 2 k -ω 2 + jω ω k Q k . ( 9
)
where k is an integer such as k ∈ [1; 2; ...; p-2; p-1], with p the total number of (null and non-null) modes. The coefficients a k , ω k and Q k are respectively the modal amplitude, the resonance pulsation and the quality factor of the k th resonance mode, and a 0 , b 0 and c 0 are the coefficients of the uniform mode.

Mathematical nature of the model

The state-of the art model is finally constituted by the three equations of the following system:

η(t) = h U j e α i W v ac (t -τ ) ∆p(t) =∆p src (t) + ∆p los (t) = ρδ d bU j W d dt tanh η(t) -y 0 b - ρ 2 v ac (t) α vc 2 sgn(v ac (t)) V ac (ω) =Y in (ω) • ∆P (ω) = a 0 b 0 jω + c 0 + p-1 k=1 a k jω ω 2 k -ω 2 + jω ω k Q k • ∆P (ω) (10)
A reformulation of system [START_REF] Gilbert | Calculation of the steady-state oscillations of a clarinet using the harmonic balance technique[END_REF] as a classical first-order system:

ẋ(t) = f (x(t), x(t -τ ), ẋ(t -τ ), λ) (11)
is helpful to highlight the mathematical nature of the model, and necessary for its implementation in DDE-Biftool. In such a formulation, x is the vector of state variables, λ the set of parameters, and f a nonlinear function. In order to improve numerical conditioning, the system under study is made dimensionless through the definition of the following dimensionless variables:

t = ω 1 t ṽ( t) = he α i W bU j v ac ( t). (12)
Finally, the components of vector x are here the variables ṽk ( t) and their first derivative with respect to time vk ( t), where ṽk ( t) is the projection of the dimensionless acoustic velocity ṽ( t) on the k th resonance mode (see equation 9).

The development is detailed in appendix A, and finally leads to the following system of 2p-1 equations (with p the total number of acoustic modes in the third equation of system [START_REF] Gilbert | Calculation of the steady-state oscillations of a clarinet using the harmonic balance technique[END_REF]):

v0 ( t) = a 0 ρδ d he α i W W b 0 p-1 i=0 zi ( t -τ ) 1 -tanh 2 p-1 i=0 ṽi ( t -τ ) - y 0 b - a 0 ρbW 2b 0 α 2 vc he α i W γ τ p-1 i=0 ṽi ( t) abs p-1 i=0 ṽi ( t) - c 0 b 0 ω 1 ṽ0 ( t) vk ( t) = zk ( t) żk ( t) = a k ρδ d he α i W W p-1 i=0 żi ( t -τ ) 1 -tanh 2 p-1 i=0 ṽi ( t -τ ) - y 0 b - 2a k ρδ d he α i W W p-1 i=0 zi ( t -τ ) 2 tanh p-1 i=0 ṽi ( t -τ ) - y 0 b 1 -tanh 2 p-1 i=0 ṽi ( t -τ ) - y 0 b - a k ρbW 2α 2 vc he α i W γ τ p-1 i=0 zi ( t) abs p-1 i=0 ṽi ( t) + p-1 i=0 ṽi ( t) p-1 i=0 zi ( t) sgn bW ω 1 he α i W γ τ p-1 i=0 ṽi ( t) - ω k ω 1 2 ṽk ( t) - ω k ω 1 Q k zk ( t) (13)
∀k ∈ [1; 2; ...; p -2; p -1].

Such a formulation highlights the presence of a delayed derivative term ż( tτ ), and thus the neutral nature of the system, already underlined in [START_REF] Auvray | Regime change and oscillation thresholds in recorderlike instruments[END_REF][START_REF] Terrien | Numerical resolution of a physical model of flute-like instruments: comparison between different approaches[END_REF].

This feature of the model constitutes an important difference with other wind instruments (such as reed instruments and brass instruments), which are modeled through ordinary differential equations.

3 Methods: orthogonal collocation and numerical continuation

Since the early works of Gilbert et al. [START_REF] Gilbert | Calculation of the steady-state oscillations of a clarinet using the harmonic balance technique[END_REF] about the determination of periodic solutions of selfsustained systems, different numerical tools have been developed or adapted to perform such calculations for musical instrument models.

Particularly, different studies have dealt with both computation and continuation of periodic solutions and bifurcation analysis of reed instruments with the software AUTO [START_REF] Karkar | Toward the systematic investigation of periodic solutions in single reed woodwind instruments[END_REF], Harmbal [START_REF] Fritz | Some aspects of the harmonic balance method applied to the clarinet[END_REF][START_REF] Farner | Contribution to harmonic balance calculations of self-sustained periodic oscillations with focus on single-reed instruments[END_REF] or Manlab [START_REF] Karkar | Oscillation threshold of a clarinet model: a numerical continuation approach[END_REF][START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]. These software differ in numerical methods they use, both for the calculation of periodic solutions (harmonic balance and collocation for Manlab, collocation for AUTO), and for the continuation algorithm (predictorcorrector method for AUTO, asymptotic numerical method for Manlab).

As underlined in the previous section and in [START_REF] Auvray | Regime change and oscillation thresholds in recorderlike instruments[END_REF], even the simplest model of flute-like instrument includes an additional delay compared to models of other wind instruments. As they do not deal with delayed differential equations, the different software mentioned above are not helpful for the present study. Among the few distributed numerical tools dealing with computation and continuation of periodic solutions in delayed systems, as Knut [START_REF] Szalai | A continuation and bifurcation software for delay-differential equations[END_REF] and DDE-Biftool [START_REF] Engelborghs | DDE Biftool: a Matlab package for bifurcation analysis of delay differential equations[END_REF], only DDE-Biftool has been subsequently adapted to treat the case of neutral delayed differential equations [START_REF] Barton | Collocation schemes for periodic solutions of neutral delay differential equations[END_REF][START_REF] Barton | Bifurcation analysis tools for neutral delay equations: a case study[END_REF]. Throughout this paper, numerical results are thus computed using DDE-Biftool and its extension for neutral systems, both based on orthogonal collocation method [START_REF] Engelborghs | Collocation methods for the computation of periodic solutions of delay differential equations[END_REF][START_REF] Barton | Collocation schemes for periodic solutions of neutral delay differential equations[END_REF][START_REF] Engelborghs | Numerical bifurcation analysis of delay differential equations using DDE-Biftool[END_REF] for computation of periodic solutions and on a predictorcorrector approach for continuation of solution branches. The general principles of orthogonal collocation, and its extension to neutral delayed equations are thus briefly recalled below.

Principles of the method

Numerical resolution of a system of ordinary, delayed or neutral differential equations implies a discrete representation of the unknown solution. Orthogonal collocation method allows to compute periodic solutions, based on a temporal discretisation of a single period. As we are interested in self-oscillating systems, it is important to note that the period T of the solution is an unknown of the problem. The principle of the method used here consists in approximating the real solution by a piecewise polynomial.

More precisely, a single (unknown) period of the periodic solution is divided in N intervals, which constitutes a mesh Π :

Π = [t 0 = 0; t 1 ; t 2 ; ...; t N -1 ; t N = T ]. ( 14
)
The number N of intervals influences on the results accuracy.

Each interval [t i ; t i+1 ] is again discretised on a set of m+1 representation points, distributed uniformly on the interval. On each interval, the solution is interpolated by Lagrange polynomials. On an interval defined as t ∈ [t i ; t i+1 ], the approximate solution x(t) of the real solution x(t) is thus defined as:

x(t) = m j=0 x i+ j m • P i,j (t) (15)
where x i+ j m is the value of the solution at the representation point t i+ j m = t i + j m (t i+1 -t i ), and P i,j (t) the Lagrange interpolating polynomial, defined as:

P i,j (t) = Π m r=0;r =j t -t i+ r m t i+ j m -t i+ r m (16)
The unknowns are thus the values of x i+ j m , and the period T . Taking into account the continuity conditions at the mesh points t i , one obtains, for a L-dimensional system, L(N m + 1) + 1 unknowns.

For each interval, the equations are not written at the m + 1 representation points, but at the m collocation points, defined as the zeros of the m-order Legendre polynomial for the considered interval. This choice of collocation points gives an optimal rate of convergence as the number of mesh points increases. Projection of the differential equations on these N × m collocation points leads to an algebraic system of N mL equations. The L + 1 additional equations required to close the system are provided by the periodicity condition: x(0) = x(T ) (L equations) and by a phase condition. Indeed, as the solutions are periodic, if x(t) is a solution, x(t+φ) is also a solution, regarless the value of φ. A phase condition is therefore introduced to fix the phase origin. A very common choice, adopted here, is the integral phase condition, initially introduced by Doedel [START_REF] Doedel | AUTO: A program for automatic bifurcation analysis of autonomous systems[END_REF]. Finally, as detailed in [START_REF] Engelborghs | Collocation methods for the computation of periodic solutions of delay differential equations[END_REF], Newton method is applied to solve the obtained algebraic system.

Adaptation of the method to neutral delay systems

The extension of orthogonal collocation method to neutral delayed differential equations imposes to define an interpolation scheme for the approximation of both the delayed and the derivative delayed variables. Various schemes have been proposed for the case of non-neutral delayed systems, as for example in [START_REF] Engelborghs | Collocation methods for the computation of periodic solutions of delay differential equations[END_REF]. For the case of neutral systems, Barton [START_REF] Barton | Collocation schemes for periodic solutions of neutral delay differential equations[END_REF] proposed the following approximations, implemented in the extension of DDE-Biftool for neutral equations:

x(t -τ ) = m j=0 x k+ j m • P k,j (t) ẋ(t -τ ) = m j=0 x k+ j m • P ′ k,j (t) (17)
where P k,j (t) corresponds to the Lagrange interpolating polynomial defined in equation ( 16), and P ′ k,j (t) to its derivative. k is an integer such that (tτ ) ∈ [t k ; t k+1 ]. For more details, we refer the reader to [START_REF] Barton | Collocation schemes for periodic solutions of neutral delay differential equations[END_REF] and [START_REF] Barton | Bifurcation analysis tools for neutral delay equations: a case study[END_REF].

Stability of static and periodic solutions

Let us consider the following system of neutral differential equations:

ẋ(t) = f (x(t), x(t -τ ), ẋ(t -τ ), λ). ( 18
)
For sake of clarity, we consider here a system with a single delay, but the following considerations can be generalised to multi-delay systems.

A static solution x s of sytem [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF] satisfies the following equation:

f (x s , x s , 0, λ s ) = 0. ( 19
)
Stability property of such an equilibrium solution is given by the stability of equation ( 18) linearised around the considered static solution x s :

ẋ(t) =A 0 [x(t) -x s ] + A 1 [x(t -τ ) -x s ] + A 2 ẋ(t -τ ), (20)
where A i is the partial derivative of the nonlinear function f with respect to its (i + 1) th argument.

Defining a new variable y(t) = x(t)-x s , equation [START_REF] Dhooge | Matcont: a matlab package for numerical bifurcation analysis of odes[END_REF] becomes:

ẏ(t) = A 0 y(t) + A 1 y(t -τ ) + A 2 ẏ(t -τ ). ( 21
)
The stability of this linearised equation depends on the value of the roots κ of its associated characteristic equation:

det κI -A 0 -A 1 e -κτ -κA 2 e -κτ = 0 ( 22
)
where I is the identity matrix. This transcendental equation has an infinite number of roots, and its resolution requires specific numerical techniques (see for example [START_REF] Breda | Pseudospectral approximation of eigenvalues of derivative operators with non-local boundary conditions[END_REF]). If all the roots κ have negative real parts, any disturbance of the equilibrium solution is attenuated with time, and the solution is locally stable. Conversely, if at least one of these eigenvalues have a positive real part, a small disturbance superimposed on the solution is amplified with time, and the solution is thus locally unstable. We precise local stability because this calculation involves the linearisation of the system around a solution, and the analysis is thus only valid for small perturbations around the considered solution. Further details and demonstrations are provided in [START_REF] Krauskopf | Numerical continuation methods for dynamical systems[END_REF][START_REF] Breda | Pseudospectral approximation of eigenvalues of derivative operators with non-local boundary conditions[END_REF].

When a stable equilibrium solution becomes unstable, the system behaviour after this bifurcation point depends on the real, pure imaginary or complex feature of the eigenvalue which crosses the imaginary axis (see for example [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF]).

Throughout this paper, we particularly focus on the case where two complex conjugate eigenvalues cross the imaginary axis. This scenario corresponds to a Hopf bifurcation, which gives rise to a periodic solution, and is thus particularly interesting in the case of musical instruments.

This local stability analysis can be generalised to periodic solutions (see [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF]), whose stability properties then depend on the eigenvalues of the monodromy operator, the so-called Floquet multipliers [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF][START_REF] Krauskopf | Numerical continuation methods for dynamical systems[END_REF]. A periodic solution is stable as long as all its Floquet multipliers lie in the unit circle. As previously, when a periodic solution loses its stability, the resulting regime observed after this bifurcation point depends on which way the Floquet multipliers leave the unit circle at the bifurcation point.

Continuation method

Starting from a static or periodic solution x 0 computed for a set of parameters λ 0 , numerical continuation methods allow (under sufficient conditions on the smoothness of the nonlinear function f ) to follow the solution when a parameter -the continuation parameter -varies. Thus, one can access to a solution family (x, λ), called branch of solutions, which contains the solution (x 0 , λ 0 ). As underlined for example in [START_REF] Karkar | Toward the systematic investigation of periodic solutions in single reed woodwind instruments[END_REF], different numerical continuation methods exist: one can cite classical predictor-corrector methods, which rely on a discrete representation of the branches, and the asymptotic numerical method [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], which attempts to compute continuous portions of the solution branches.

The software DDE-Biftool, used for the present study, uses a predictor-corrector continuation algorithm: knowing a point (x 0 , λ 0 ) of a branch, the neighboring point (for a set of parameter λ 0 + ∆λ) is first estimated through a secant predictor [START_REF] Krauskopf | Numerical continuation methods for dynamical systems[END_REF][START_REF] Engelborghs | DDE Biftool: a Matlab package for bifurcation analysis of delay differential equations[END_REF], and then corrected using an iterative Newton-Raphson algorithm. This process thus allows to compute the whole solution branch.

To improve numerical conditioning, and avoid the common but pathological case of a vertical tangent -which corresponds to a fold (or saddlenode) bifurcation [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF] -the solution branch is parameterised through the Keller pseudo-arclength equation, which can be seen as a kind of curvilinear abscissa. For more details, we refer the reader to [START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF].

Numerical continuation algorithm, coupled to the orthogonal collocation method, can be used to calculate bifurcation diagrams of the system of interest. In an ideal case, such diagrams show all the static and periodic solutions as a function of the continuation parameter. Through knowledge of unstable solutions, coexistence of solutions, and bifurcations of periodic solutions, it provides a more global vision of the system dynamics than time-domain simulations alone. However, computation of a solution branch through numerical continuation requires the knowledge of at least two points of this branch (i.e. two starting points) to enable the secant prediction of an other point, which raises the question of initialisation.

For the model studied here, the static solution x s = 0 being known analytically, the initialisation of the corresponding branch is not problematic. The interest of numerical continuation then lies in the computation of the eigenvalues κ (see section 3.3) along this branch, and thus in the detection of the Hopf bifurcations.

At each bifurcation point, one can follow the emerging periodic solution branch. One of the two required starting points is constituted by the degenerate periodic solution of zero amplitude, which corresponds to the Hopf point itself. The second starting point is first predicted as a monochromatic periodic solution of small amplitude, whose frequency corresponds to that given by the imaginary part of the critical eigenvalue at the Hopf point. This solution is then corrected through the Newton-Raphson algorithm. However, this method does not ensure the convergence of the correction process. A possible improvement should rely on a new way of prediction of the starting solution, using the eigenvector related to the critical eigenvalue, as done in DDE-Biftool for non neutral systems [START_REF] Engelborghs | Numerical bifurcation analysis of delay differential equations using DDE-Biftool[END_REF].

Due especially to the fact that we have only access to an approximation of the eigenvalues κ (and thus to the different Hopf bifurcations), and to the possible existence of isolated periodic solutions (i.e. not connected to the static solution branch), it remains impossible to be sure that a complete bifurcation diagram has been obtained. [START_REF] Schumacher | Self-sustained oscillations of organ flue pipes: an integral equation solution[END_REF] Application to flute-like instruments

Comparison of the results of collocation and continuation with timedomain simulations

In order to validate the continuation approach, we first focus on the comparison, in terms of amplitude, frequency and waveforms, between the results of a classical time-domain solver based on a Runge-Kutta method on the one hand and the bifurcation diagram on the other hand.

As it is directly related to the pressure in the musician's mouth (and thus to the choice of whether that he blows hard or not), the dimensionless delay τ is a particularly interesting parameter, and is chosen as bifurcation parameter.

Paramater values used for the different calculations are summarised in table 1. We consider here the ideal case of a cylindrical resonator without side holes, 400 mm in length and 16 mm diameter, which is close to the typical dimensions of an alto recorder. The contribution of the excitation window is taken into account through a correction term in the analytical formula of the input admittance (see [START_REF] Chaigne | Acoustique des instruments de musique (Acoustics of musical instruments)[END_REF] for more details). From the input admittance calculated analytically, a fitting algorithm allows the extraction of the modal coefficients of the six first resonance modes (a uniform mode at zero frequency, and 5 non uniform modes). The comparison, in terms of modulus and phase, between the original and the fitted admittances, is provided in figure 5. Between these two admittances, the relative deviations in terms of amplitude of the resonance peaks are 12, 2% for the fifth non uniform mode, 2, 1% for the fourth non uniform mode, and lower than 1% for the three first non uniform modes. In terms of resonance frequencies, the relative deviations between the original and the fitted admittances are lower than 0, 1% for all the five non uniform modes.

One can note that an increase of the number of modes taken into account allows a better modeling of the input admittance. However, the size of the system being directly related to the number of modes, adding modes considerabily increases the computation cost. Moreover, on such instruments, most of the notes are played on the two first registers, which corresponds to oscillations on the two first resonance modes. Since the two first (non uniform) modes are correctly fitted by taking into account six resonance modes (see figure 5), this number seems to be a reasonable compromise.

For sake of clarity, we first consider a single periodic solution branch, corresponding to the second register, that is to say to oscillations at a frequency close to the second resonance frequency of the input admittance.

Figure 6 represents the amplitude of the acoustic velocity v ac as a function of the dimensionless delay τ , for both the results obtained through time-domain integration for three fixed τ values (τ = 0.6, τ = 0.846 and τ = 1.2) and the branch computed thanks to numerical continuation algorithm. Since it leads to convergence of the solver and permits an easy resampling of the signals at a sampling frequency suitable for audio production systems, a sampling frequency f s = 10 × 44100 Hz is used for time-domain simulations. Orthogonal collocation is achieved using 75 mesh points per period, and an interpolation polynomial of degree 5. This comparison shows good agreement between the results of the two methods, with a maximum relative deviation of the amplitude of Figure 5: Analytically calculated and fitted input admittances of a cylindrical resonator of flutelike instrument, 400 mm in length and 16mm diameter. Relative deviations between the curves in terms of amplitude of the resonance peaks: 12.2% for the fifth non uniform mode, 2.1% for the fourth, and lower than 1% for the other modes. In terms of resonance frequencies, the relative deviations between the two curves are lower than 0.1 % for all the modes taken into account. 0.36 % at τ = 1.2.

Figure 7 provides, for the same data, the comparison in terms of oscillation frequency. As for the amplitude, this comparison shows good agreement, with a maximum relative deviation of 0.1 % at τ = 0.6, which corresponds to 1.7 cent.

In the same way, figure 8 provides the comparison between the two methods, for both the calculated waveforms and amplitude of the 15 first harmonics of the spectrum, at the three different points of the branch already represented in figure 6. As the fitted input admittance used here takes into account six resonance modes, one can wonder about the physical meaning of the higher harmonics represented in this figure. However, the aim of this comparison is here to check that the two methods converge to a same solution. In the case of a study focusing on the sound spectrum, one should probably consider a larger number of resonance modes. Although the sampling frequency used for simulations and the number of discretisation points retained for orthogonal collocation are not "equivalent" (and thus the two methods do not resolve the same number of harmonics), the comparisons provided in figure 8 show again good agreement. Especially, the calculation of the L2-norm of the different waveforms highlights rel-

0 (m 2 • kg -1 ) 642 b 0 83 c 0 (s -1 ) 294 a 1 (m 2 • kg -1 ) 15 ω 1 (rad • s -1 ) 2480 Q 1 50 a 2 (m 2 • kg -1 ) 14 ω 2 (rad • s -1 ) 4982 Q 2 64 a 3 (m 2 • kg -1 ) 12 ω 3 (rad • s -1 ) 7501 Q 3 72 a 4 (m 2 • kg -1 ) 10 ω 4 (rad • s -1 ) 10040 Q 4 77 a 5 (m 2 • kg -1 ) 9 ω 5 (rad • s -1 ) 12595 Q 5 77 α i (m -1 ) 0.4/h = 400 δ d (m) 4 π (2hW ) = 0.0036 b (m) 2h/5 = 0.0004 y 0 (m) 0.0001 h (m) 0.001 W (m) 0.00425 ρ (kg • m -3 ) 1.2 c p (m • s -1 ) 0.4U j α vc 0.6
Table 1: Parameter values used for numerical resolutions of the model.

ative deviations between the results of orthogonal collocation and time-domain simulation of 0.78% at τ = 0.6, 0.46% at τ = 0.846 and 0.04% at τ = 1.2. Although these results focus on a single solution branch, and thus lead to a very partial knowledge of the model behaviour, they nevertheless provide, compared to more classical analysis and resolution methods (such as linear analysis and time-domain simulations), new information about the functioning of the instrument. For sake of consistency with the literature, oscillation amplitude and frequency (already represented respectively in figures 6 and 7) are represented in figures 9 and 10 as functions of the reduced jet velocity θ: Such representations shed some light on different known behaviours of the instrument. Especially, it highlights that the methods used here predict precisely the saturation of the oscillation amplitude (figure 9), a commonly observed behaviour in experiments and simulations [START_REF] J W Coltman | Time-domain simulation of the flute[END_REF][START_REF] Auvray | Regime change and oscillation thresholds in recorderlike instruments[END_REF][START_REF] Fletcher | Sound production by organ flue pipes[END_REF][START_REF] Meissner | Aerodynamically excited acoustic oscillations in cavity resonator exposed to an air jet[END_REF], which is not explained or predicted by the often-used linear analysis of the model [START_REF] Auvray | Regime change and oscillation thresholds in recorderlike instruments[END_REF]. In the same way, if the strong dependance of the frequency on the jet velocity U j , highlighted in figure 10, is a well-known behaviour of both models and real instruments (see for example [START_REF] J W Coltman | Time-domain simulation of the flute[END_REF][START_REF] Auvray | Regime change and oscillation thresholds in recorderlike instruments[END_REF][START_REF] Terrien | Flutelike musical instruments: A toy model investigated through numerical continuation[END_REF][START_REF] Meissner | Aerodynamically excited acoustic oscillations in cavity resonator exposed to an air jet[END_REF][START_REF] J W Coltman | Resonance and sounding frequencies of the flute[END_REF][START_REF] Giordano | Direct numerical simulation of a recorder[END_REF]), a linear analysis of the model only gives a rough estimation of the frequency evolution, and does not distinguinsh between stable and unstable parts of the branch. As highlighted in figure 10, the bifurcation diagram not only predicts precisely the frequency evolution along the branch, but also the stabilisation of the frequency slightly above the resonance frequency (observed experimentally for example in [START_REF] Auvray | Regime change and oscillation thresholds in recorderlike instruments[END_REF]). Through the computation of stability properties of the branch of periodic solutions, it finally gives information about the minimum and maximum frequencies that can be observed for a given periodic regime. Thereby, figure 10 highlights that the oscillation frequency of the second register can evolve between 776 Hz (at θ = 12) and 804 Hz (at θ = 34). By contrast, as a linear analysis only provides information about the instability of the static solution, which corresponds to the the birth of the periodic solution branch, it would merely gives an estimation of the frequency at the branch extremities (here at θ = 10.5 and θ = 152).

θ = U j • 2π W ω 1 . ( 23
)
To conclude, this first study of a single branch of periodic solutions allows to validate, for the state-of-the art model of flute-like instruments, the approach through orthogonal collocation and numerical continuation, and to begin to appreciate the benefits of such an approach. Hereafter, we extend the study to a more complete bifurcation diagram, in order to provide a better understanding of the instrument behaviour.

Analysis of the transition between registers

The phenomenon of "register change" is wellknown by flute players: occuring while the musician blows harder in the instrument, it corresponds to a "jump" from a given note on the first register to another one on the second register whose frequency is approximately twice the first (that is to say, a jump of approximately one octave).

As highlighted in figure 11, representing the oscillation frequency of the acoustic velocity v ac (t) as a function of the dimensionless delay τ , this phenomenon is also observed with time-domain simulations of the model, where the parameter τ now varies linearly. Starting from a frequency of 400 Hz at τ = 1.2, a decrease of the delay (corresponding to an increase of the jet velocity U j ) leads to successive jumps to notes with higher frequencies. As observed experimentally (see for example [START_REF] Auvray | Regime change and oscillation thresholds in recorderlike instruments[END_REF]), the register changes are accompanied by hysteresis: starting from the arrival point of the decreasing ramp of delay (at τ = 0.1) and performing an increasing ramp of the delay also leads to jumps, but they occur at different delay values from those observed previously.

However, if this approach highlights the presence of register changes with hysteresis, it does not provide information about the cause of these jumps and of the associated hysteresis phenomenon. Moreover, as can be seen in figure 12, representing for the same simulations as in figure 11 the oscillation amplitude as a function of the dimensionless delay, the simulation results can be difficult to understand and interpret, due to the large number of observed regimes and "amplitude jumps".

The sensitivity to initial conditions make such an interpretation even more complicated: for the same parameters, a small change in the initial conditions of the state variables can result in a very different behaviour (periodic solution with different amplitude and frequency, but also non- Oscillation frequency obtained through time-domain integration of the flute model for both a decreasing and an increasing ramp of the delay τ . The lines of blue dots indicate the resonance frequencies of the resonator The superimposition, in figures 13 and 14 of the bifurcation diagram and the results of simulations already presented above, provides a global overview of the different existing periodic regimes, and thus allows a better understanding of the observed phenomena. The bifurcation diagram highlights, for the range of the delay τ under investigation, the existence of five periodic solution branches, corresponding to the five different registers, that is to say to oscillations on each (non uniform) mode of the resonator (see the third equation of system [START_REF] Gilbert | Calculation of the steady-state oscillations of a clarinet using the harmonic balance technique[END_REF]). Stable and unstable parts of the branches are respectively represented by solid lines and dashed lines. Oscillation frequency obtained through time-domain integration of the flute model for both a decreasing and an increasing ramp of the delay τ , superimposed with the bifurcation diagram obtained through numerical continuation.

(f 1 = ω 1 2π , f 2 = ω 2 2π , f 3 = ω 3 2π , f 4 = ω 4 2π , f 5 = ω 5 2π ).
This comparison between the bifurcation diagram and simulations shows that the "jump" from a given register A to another register B is caused by the loss of stability (through a Neimark-Sacker bifurcation) of the branch of periodic solution corresponding to the register A. Such a bifurcation can lead to the birth of a quasiperiodic regime, which would be called a multiphonic sound in a musical context, and whose generation mechanism in flute-like instruments is described more precisely in [START_REF] Terrien | Flutelike musical instruments: A toy model investigated through numerical continuation[END_REF].

The confrontation between the bifurcation diagram and simulations also allows to explain the hysteresis phenomenon by the coexistence of several stable periodic solutions for a same value of the delay τ . Let us consider the case of transition between registers 4 and 5 (corresponding respectively to oscillations at about ω 4 /(2π) = 1600Hz Oscillation amplitude obtained through time-domain simulation of the flute model for both a decreasing and an increasing ramp of the delay τ , superimposed with the bifurcation diagram obtained through numerical continuation.

and ω 5 /(2π) = 2000Hz): starting from the 4th register, at τ = 0.44, a decrease of τ leads to a jump on the five register at τ = 0.19. Starting an increasing ramp of the delay from this point, one observe that the system remains on the branch corresponding to the 5th register, until it becomes unstable, at τ = 0.45. In the same way, it can explain the fact that some regimes (for example the third register) are only observed, using timedomain integration, in the case of an increasing ramp of the delay, and not in the case of a decreasing ramp (see figures 11 and 13).

Thereby, the bifurcation diagram allows us to predict, in the case of a quasi-static variation of the delay, both oscillation threshold and thresholds of register change. As such a prediction is impossible with the classical approach involving a linear analysis of the model (see for example [START_REF] Auvray | Regime change and oscillation thresholds in recorderlike instruments[END_REF]), it highlights the valuable contribution of the method.

As an illustration, figures 15 and 16 respectively represent, as functions of the reduced jet velocity θ, the same bifurcation diagrams as in figures 13 and 14, superimposed with the thresholds of the different periodic solution branches predicted through an analysis of the model linearised around its equilibrium solution. These thresholds correspond to the Hopf bifurcation points, whose detection is described in section 3.3. As highlighted in figures 15 and 16, such a linear analysis allows to predict the thresholds corresponding to instabilities of the equilibrium solution (that is to say the birth of the periodic solution branches), but it does not provide any information about stability of the emerging periodic solution. As a periodic solution is not necessary stable at the birth of the branch, the threshold predicted through a linear analysis does not necessary correspond to the "observable" threshold of a given oscillation regime. As a consequence of thresholds prediction and stability analysis, a bifurcation diagram also predicts and explains hysteresis ranges, which is impossible through a linear analysis, and not always possible using time-domain simulations, due to the sensitivity to initial conditions, as underlined above. This is however a relevant information in the context of musical acoustics since it is related with how confortable a musician feels when playing the instrument [START_REF] Auvray | Regime change and oscillation thresholds in recorderlike instruments[END_REF]. Indeed, an important hysteresis allows him to play both forte on the first register and piano on the second register. W ω 1 , and superimposed with the thresholds predicted through a linear analysis of the model. The additional window (a) represents a zoom on the branch related to the second register: it highlights the difference between the thresholds predicted through linear analysis, and the "observable" thresholds, corresponding to the beginning of the stable part of the branch (in solid line).

Thus, the bifurcation diagram, giving access to new kind of information (especially existence of unstable solutions, co-existence of multiple stable solutions, and oscillation amplitude and frequency of different periodic solutions of the nonlinear model), enables an easier interpretation of the simulation results (see especially figures 12 and 14), a more precise prediction of some char- W ω 1 , and superimposed with the thresholds predicted through a linear analysis of the model. The additional window (a) represents a zoom on the branch related to the second register: it highlights the difference between the thresholds predicted through linear analysis, and the "observable" thresholds, corresponding to the beginning of the stable part of the branch (in solid line). acteristics of the different periodic regimes (especially their thresholds), and explains some aspects of the model behaviour observed in simulations, such as hysteresis associated to register changes.

Conclusion

Widely used in many scientific fields, methods dedicated to computation of periodic solutions of nonlinear dynamical systems have demonstrated, in several previous works [START_REF] Gilbert | Calculation of the steady-state oscillations of a clarinet using the harmonic balance technique[END_REF][START_REF] Fritz | Some aspects of the harmonic balance method applied to the clarinet[END_REF][START_REF] Farner | Contribution to harmonic balance calculations of self-sustained periodic oscillations with focus on single-reed instruments[END_REF][START_REF] Karkar | Toward the systematic investigation of periodic solutions in single reed woodwind instruments[END_REF][START_REF] Karkar | Oscillation threshold of a clarinet model: a numerical continuation approach[END_REF], their benefits for the study of musical instruments. However, no available software (such as AUTO [START_REF] Doedel | AUTO: A program for automatic bifurcation analysis of autonomous systems[END_REF], Manlab [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], or DDE-Biftool [START_REF] Engelborghs | DDE Biftool: a Matlab package for bifurcation analysis of delay differential equations[END_REF]) is suitable for neutral delayed systems, involved in the modeling of flute-like instruments. The orthogonal collocation method has recently been adapted to such systems [START_REF] Barton | Collocation schemes for periodic solutions of neutral delay differential equations[END_REF], allowing us to consider a study of the state-of-the-art model of flutes through numerical continuation.

The comparison between the results of orthogonal collocation coupled with numerical continuation algorithm and time-domain integration shows very good agreement in terms of oscillation amplitude and frequency, thresholds of the different periodic regimes, and waveforms. It thus allows the validation of this new approach to study flute-like instruments. Giving access to bifurcation diagram of the model, this method provides a more global knowledge of the model dynamics, and permits to explore more broadly the influence of the variation of a parameter. Thereby, it allows to explain and predict more precisely different aspects of the model behaviour, which can be related to some phenomena observed experimentally in flute-like instruments, such as register changes and associated hysteresis. Thus, the access to this more global knowledge of the model behaviour should allow to envisage an easier qualitative confrontation between numerical and experimental results.

Nevertheless, this method presents a number of restrictions, and is thus not capable of replacing time-domain simulations entirely, but is a useful complement. Especially, one can only compute steady-state solutions, and there is no calculation of the transients, which are particularly important in the musician's playing. In the same way, this method only provides information about the behaviour of the model in the case of a quasi-static variation of the continuation parameter. Since recent studies [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] have highlighted the strong influence of the variation rate of a control parameter on the behaviour of a dynamical system, one can wonder about the validity of the bifurcation diagram in a musical context, where many control parameters vary continuously. Moreover, while orthogonal collocation can be extented to compute non-periodic solutions (such as for example quasiperiodic, intermittent or chaotic regimes), the software package used here is not currently suitable for such calculations.

A Reformulation of the model

Both to highlight its mathematical nature and to allow its implementation in DDE-Biftool and its extension for neutral systems, the model needs to be written as a first-order system: ẋ(t) = f (x(t), x(tτ ), ẋ(tτ ), λ), with x the vector of state variables, and λ the set of parameters.

Starting from system [START_REF] Gilbert | Calculation of the steady-state oscillations of a clarinet using the harmonic balance technique[END_REF], and writting the third equation in the time-domain leads to:

a k jω ω 2 k -ω 2 + jω ω k Q k • ∆P (ω) (24)
For each non uniform mode k, we define V ac k (ω) such as V ac k (ω) =

a k jω ω 2 k -ω 2 +jω ω k Q k
• ∆P (ω)(with k an integer such as k ∈ [1; p -1]). In the same way, we define V ac 0 (ω) = a 0 b 0 jω+c 0 • ∆P (ω), such as one finally obtains V ac = V ac 0 + p-1 k=1 V ac k . System [START_REF] Chaigne | Acoustique des instruments de musique (Acoustics of musical instruments)[END_REF] can thus be written as: V ac k (ω)

η(t) = h U j e α i W v ac (t
V ac 0 (ω) = a 0 b 0 jω + c 0

• ∆P (ω)

V ac 1 (ω) = a 1 jω ω 2 1 -ω 2 + jω ω 1 Q 1 • ∆P (ω)
. . .

V ac (p-1) (ω) =   a (p-1) jω ω 2 (p-1) -ω 2 + jω ω (p-1) Q (p-1)   • ∆P (ω) (25)
Through an inverse Fourier transform, each expression V ac k (with k ∈ [1; p -1]) can be written as:

d 2 v ac k (t) dt 2 + ω k Q k • dv ac k (t) dt + ω 2 k v ac k (t)
= a k d∆p(t) dt [START_REF] De | The sound of oscillating air jets: Physics, modeling and simulation in flute-like instruments[END_REF] and the expression of V ac 0 becomes: b 0 • dv ac 0 (t) dt + c 0 v ac 0 (t) = a 0 ∆p(t) [START_REF] Rayleigh | The theory of sound second edition[END_REF] The reinjection of the expressions of ∆p(t) and η(t) leads to:

v ac (t) = p-
- ω k Q k dv ac k (t) dt -ω 2 k v ac k (t). (28)
To improve numerical conditioning of the problem, we define the following dimensionless variables:

t =ω 1 t ṽk ( t) = he α i W bU j v ac k ( t) (29)
Noting v( t) and v( t) respectively the first and second order derivative of ṽ( t) with respect to the dimensionless time t, one finally obtains:

- ω k ω 1 Q k vk ( t) (30)
∀k ∈ [1; 2; ...; p -2; p -1].

The definition of the variables: zk ( t) = vk ( t) finally allow to reduce system [START_REF] J W Coltman | Jet drive mechanisms in edge tones and organ pipes[END_REF] to a first-order system:

- ω k ω 1 Q k zk ( t) (31)
∀k ∈ [1; 2; ...; p -2; p -1], with p the total number of modes.

Figure 1 :

 1 Figure 1: Schematic representation of the principle of numerical continuation. x is a given characteristic of the periodic solution, as for example its amplitude or its frequency. (a): periodic solutions obtained through, for example, orthogonal collocation, represented as a fonction of the parameter λ. (b): branch of periodic solution obtained through numerical continuation.

Figure 2 :

 2 Figure 2: Recorder section and simplified representation of its exciter, constituted by the interaction between an unstable air jet created at the channel exit, and a sharp edge called labium.

Figure 3 :

 3 Figure 3: Schematic representation of the jet behaviour, based on Fabre in [24]. (a) Perturbation of the jet at the channel exit by the acoustic field present in the resonator. (b) Convection and amplification of the perturbation, due to the unstable nature of the jet. (c) Jet-labium interaction: oscillation of the jet around the labium, which sustains the acoustic field.

Figure 6 :

 6 Figure6: Oscillation amplitude as a function of the dimensionless delay τ for a single periodic solution branch : comparison between the results of numerical continuation and time-domain integration. The maximum relative deviation between the two methods is 0.36% at τ = 1.2.

2 Figure 7 :

 27 Figure7: Oscillation frequency as a function of the dimensionless delay τ for a single periodic solution branch : comparison between the results of numerical continuation and time-domain integration. The maximum relative deviation between the two methods is 0.1% at τ = 0.6. The dark dot-dashed line indicates the second resonance frequency (f 2 = ω 2 2π ) of the resonator.

Figure 8 :

 8 Figure 8: Comparison between the results of orthogonal collocation and time-domain integration in terms of waveform (left) and spectrum magnitude (right) represented in dB (reference amplitude: V ref = 1). Each line corresponds to a value of τ (τ = 0.6 for the first line, τ = 0.846 for the second line, τ = 1.2 for the third line). In terms of L2-norm of the different waveforms, relative deviations between the results of the two methods are of 0.78% at τ = 0.6, 0.46% at τ = 0.846 and 0.04% at τ = 1.2.

Figure 9 :

 9 Figure 9: Oscillation amplitude as a function of the reduced jet velocity θ for a single periodic solution branch. Superimposition of the results of numerical continuation and time-domain integration.

2 Figure 10 :

 210 Figure 10: Oscillation frequency as a function of the reduced jet velocity θ for a single periodic solution branch. Superimposition of the results of numerical continuation and time-domain integration. The dark dot-dashed line indicates the second resonance frequency (f 2 = ω 2 2π ) of the resonator.

  Figure 11:Oscillation frequency obtained through time-domain integration of the flute model for both a decreasing and an increasing ramp of the delay τ . The lines of blue dots indicate the resonance frequencies of the resonator(f 1 = ω 1 2π , f 2 = ω 2 2π , f 3 = ω 3 2π , f 4 = ω 4 2π , f 5 = ω 5 2π).

  limit cycle of Vac (peak to peak) simulation -decreasing ramp of τ simulation -increasing ramp of τ

Figure 12 :

 12 Figure 12:Oscillation amplitude obtained through time-domain integration of the flute model for both a decreasing and an increasing ramp of the delay τ .

  Vac (Hz) continuation -stable solutions continuation -unstable solutions simulation -decreasing ramp of τ simulation -increasing ramp of τ

Figure 13 :

 13 Figure 13:Oscillation frequency obtained through time-domain integration of the flute model for both a decreasing and an increasing ramp of the delay τ , superimposed with the bifurcation diagram obtained through numerical continuation.

  Figure 14:Oscillation amplitude obtained through time-domain simulation of the flute model for both a decreasing and an increasing ramp of the delay τ , superimposed with the bifurcation diagram obtained through numerical continuation.

Figure 15 :

 15 Figure 15: Oscillation frequency obtained through numerical continuation, represented as a function of the reduced jet velocity θ = U j •2π

Figure 16 :

 16 Figure 16: Oscillation amplitude obtained through numerical continuation, represented as a function of the reduced jet velocity θ = U j •2π

  v0 ( t) = a 0 ρδ d he α i W a k ρδ d he α i W W

-

  v0 ( t) = a 0 ρδ d he α i W =z k ( t) żk ( t) = a k ρδ d he α i W W 2a k ρδ d he α i W W


Emmanuel Hugonnot

Marie-Hélène Delville

Jean-Pierre Delville
email: jp.delville@cpmoh.u-bordeaux1.fr

Dissymmetrization of Microparticle Surface by Laser-Induced Photochemical Deposition

Keywords: numbers: 81.15.Fg, 81.10.Dn, 42.40.Eq, 42.70.Gi 2 photochemical deposition, surface patterning, nano/micro-sized particles, laser, chromate

come

Dissymmetrization of micro-particle surface by laser-induced photochemical deposition

Introduction

Nano-and micro-size particles have an increasingly important role in emerging nano/micro "technologies". Many efforts have been made to prepare isotropically coated and functional [1] colloid particles. Surprisingly, much less attention has been dedicated to the dissymmetrization of their chemical properties, since dissymmetry brings an additional degree of freedom that can be used to increase the complexity of targeted micro/nanomaterials. Such a strategy could be very promising for self-assembling, [START_REF] Lehn | Supramolecular Chemistry: Concepts and Perspectives[END_REF], chemical sensoring, [START_REF] Himmelhaus | [END_REF], drug delivery, [4], or micro-optics [5]. However, conventional and advanced lithographic techniques [6] are essentially limited to planar and cylindrical substrates [7]. The first attempts for particles dissymmetrisation were investigations of partial silanization [8], chemisorption in directional fluxes [9] and differential nano-clusters adsorption at interfaces [10]. We extend these methods by presenting a general and convenient alternative that (i) combines micro-photochemical deposition driven by a strongly localized laser excitation, to "nucleate" a deposit, i.e. a dissymmetry, on the surface of a micro-sphere and (ii) control its growth by irradiation parameters. When properly rescaled, the coating growth laws measured for different excitations reduces into a single master curve helping to control and predict the desired dissymmetry. The versatility of the technique and its ability to tailor micro-scale patterning on only one of the two hemispheres are also described.

Experimental Section

Laser chemical processing is a powerful non-contact tool for micro-machining and micro-processing of materials [START_REF] Bäuerle | Laser Processing and Chemistry[END_REF]. Photo-chemically stimulated deposition represents nowadays an attractive alternative for direct material writing without using mask-based photolithography. The technique uses UV or visible photons to efficiently break chemical bonds into molecules and induce rapid thin film deposition via very moderate beam intensities. Furthermore, laser light offers localized processing by accurate spatial excitation and spectral selectivity in molecule activation. When performed in liquids, photochemical deposition can be applied to a broad range of precursors, including molecular compounds.

Moreover, since the light energy conversion yields are much higher than in gas phase, one requires relatively low power lasers for processing [START_REF] Peled | [END_REF].

Based on the photo-reduction of ( )

Cr VI ions into ( )
Cr III , used in dichromated gelatins for optical recording [ 13 ], the experiment, shown in Figure 1, is based on the photochemical deposition of a chromium oxide layer [ 14 ]. It is performed at room

Cr V ( ) Cr V ( )
Cr III [13], see Eq. 1:

ν + * Cr(VI ) h Cr(VI ) , (1a)
+ → → → * " "

Cr(VI ) ROH Cr(V ) ... Cr( III ).

(1b)

Results

According to

Discussion

To interpret these data, we use a droplet growth model [15] considering the case of small normalized optical absorption σ << , the radial growth rate of the deposit is given by:

( ) ( ) ρ ρ ρ τ = - ⋅ 2 2 0 d exp 2 I 2 , d (2)
where ( ) ∝ is predicted from Eq. 2. Thus, our approach clearly suggests a unified behavior of the dynamics [16].

The method can easily be extended to micro-scale patterning of particles by scanning the cell holder stages [17]. One mode consists in discontinuously realizing micro-patterns dot by dot under a finite exposure for each deposit followed by a substrate translation perpendicular to the beam axis (Figure 3a). In the second mode, the laser writes continuous lines instead of dots (Figure 3b).

Applications of asymmetric patterning of micro-spheres can be applied in many emerging areas. A direct application is to tailor the phase and the amplitude of micro-lens coatings in optoelectronics devices [18]. Another feature is the chemical anisotropy resulting from dissymmetrization of the micro-particle hemispheres. Indeed, the induced coating can be chemically modified with organic functions or used as a mask to modify the uncoated surface.

In a second step, the deposit can be removed, allowing a further different functionalization of the previously covered surface [10]. Non-symmetrical micro-spheres can also be used to force self-assembling in a preferred way and build organized architectures, such as photonic crystals.

Conclusion

We used laser-assisted micro-photochemical deposition to coat and dynamically tailor the surface dissymmetry of micrometric silica beads. Considering the agreement observed between the experiments and predictions, as well as the description of the deposit growth rate in scaling terms, our method can easily be extended to the photochemical deposition of noble metals [19] or semi-conductors [20] on various particles such as ,

2 SiO 2 3
Al O and . We believe that micro-photodeposition can be considered as an interesting novel method to be used for micro-lens dissymmetrization, microsphere functionalization and micrometric selfassembly.

Figures Caption:

temperature on 10

 10 NUCLEOSIL µm ® silica beads dispersed in a liquid mixture composed of potassium chromate, ethanol, hydrochloric acid and water. The samples are enclosed in a tight very thin home-made cell, composed of two glass slides separated by 30 thickness mylar spacers, to prevent disturbing thermal effects induced by the light absorption, such as convection or liquid boiling. The upper slide is silanized to eliminate any photochemical deposition on the entrance face of the cell. The cell is horizontally mounted to keep the beads at rest on the bottom plate. A computer controls translation stages, with a step accuracy of 0.5 µm, to move the cell holder in the ]. At = pH 3 , is the major species in solution. As this hexavalent form of chromium can be activated by light in the blue-green wavelength range, we used a Ar + laser ( nm 514 = λ ). The addition of an organic compound, such as an alcohol ( ), enhances the photo-redox process of ROH 4 HCrO -to amorphous chromium hydroxide ( ) 3 Cr OH . This ( ) Cr III species has a very low solubility. The major pathway from ( ) Cr VI to ( ) Cr III starts with a light-induced excitation of ( ) Cr VI followed by reduction to a intermediate and a dark reaction from to ( )

  Eqs. (1a-b), the ( ) Cr III concentration in solution increases with laser intensity. The nucleation of the precipitate occurs as soon as the solubility limit is achieved. Consequently, a photo-deposit of circular shape starts to grow onto the silica particle hemisphere intercepted by the exciting beam. See the temporal sequence of spot growth for in Figure 1. The side view of the central bead illustrated in the experimental setup clearly shows that the deposited dissymmetry is preserved until the complete covering of one hemisphere. Considering the homogeneity of the induced deposit, we performed measurements of the photodeposit growth rate. As the substrate is curved, and growth mechanisms involve solute flux along the optical direction, the image processing approximates the measured cross section of the covered hemisphere by a circular area of radius ( t instead of considering the hemispherical surface covered by the deposit. The inset of Figure 2 shows measurements, carried out with a beam waist and incident laser beam power in the range µW . At a given time, the growth rate observed increases for increasing beam power.

0 a 1 σ

 1 << , where σ is the optical absorption of the chromate solution. For 0 a 1

for 1 ρ

 1 for the reaction and ϑ is the contact angle of the deposit on the substrate. Consequently, the entire dependence in exciting beam parameters is transferred to both the rescaled deposit radius and the new thermodynamic time scale. The apparently different growth laws obtained for various conditions of optical excitations (Inset of Figure2) point out for a single-scaled dynamic when the radius of the deposit is This data reduction on a master curve is shown in Figure2.Note that a simple asymptotic regime emerges for >> , we find ρ τ ∝ , in good agreement with the observations as a power law fit of the data leads to . This results shows that the growth simply depends on the energy deposited in the medium when the deposit radius becomes much larger than that the excited zone. Agreement is also observed at smaller

Figure 1 :Figure 2 : 1 Figure 3 :

 1213 Figure 1: Schematics of the experimental setup used for investigating the radial growth of a

  Figure 3

HCrO -

TiO