Minimal blowing pressure allowing periodic oscillations in a simplified reed musical instrument model: Bouasse-Benade prescription assessed through numerical continuation

A reed instrument model with N acoustical modes can be described as a 2N dimensional autonomous nonlinear dynamical system. Here, a simplified model of a reed-like instrument having two quasi-harmonic resonances, represented by a four dimensional dynamical system, is studied using the continuation and bifurcation software AUTO. Bifurcation diagrams of equilibria and periodic solutions are explored with respect to the blowing mouth pressure, with focus on amplitude and frequency evolutions along the different solution branches. Equilibria and periodic regimes are connected through Hopf bifurcations, which are found to be direct or inverse depending on the physical parameters values. Emerging periodic regimes mainly supported by either the first acoustic resonance (first register) or the second acoustic resonance (second register) are successfully identified by the model. An additional periodic branch is also found to emerge from the branch of the second register through a period-doubling bifurcation. The evolution of the oscillation frequency along each branch of the periodic regimes is also predicted by the continuation method. Stability along each branch is computed as well. Some of the results are interpreted in terms of the ease of playing of the reed instrument. The effect of the inharmonicity between the first two impedance peaks is observed both when the amplitude of the first is greater than the second, as well as the inverse case. In both cases, the blowing pressure that results in periodic oscillations has a lowest value when the two resonances are harmonic, a theoretical illustration of the Bouasse-Benade prescription.


Introduction
An important goal of the acoustics of wind instruments is to understand key components of intonation and also the ease of playing. From the physics modelling point of view, it is interesting to study the main variables that control the playing frequency (for intonation) and the minimum mouth pressure to achieve auto-oscillations (for ease of playing). It is assumed that part of the musician's judgement of ease of playing of a note is inversely related to the sounding resistance represented by the threshold blowing pressure. Support for this hypothesis is offered by measurements on the performing properties of saxophone reeds by [1]: a significant correlation was found between the soft-hard scale on which the sounding resistance of different reeds was judged by saxophonists and the threshold pressures measured in the mouths of the performers. A useful overview of the acoustics of reed and lip wind instruments can be found in books such as [2][3][4][5][6].
It is often commented that the flaring bore of brass instruments are designed such that the input impedance are, as close as possible, harmonically related. While this alignment is said to be important for intonation, it is also likely to determine the oscillation threshold and therefore improve the ease of playing. Here, the necessity of an alignment in a harmonic series is called the "Bouasse-Benade prescription" because of what Benade wrote in his famous book [2], or in [7]: The usefulness of the harmonically related air column resonances in fostering stable oscillations sustained by a reed-valve was first pointed out by the French physicist Henri Bouasse in his book "Instruments à vents" [8]. In order to illustrate this prescription, a horn was designed to provide an air column whose resonance frequencies (frequencies of maximum input impedance) were chosen to avoid all possible integer relations between them, called "tacet horn" in [7]. The purpose of this instrument is to deliberately make the conditions for oscillation unfavorable.
The effect of inharmonicity of the two first resonance frequencies on both tone colour and ease of playing have been examined experimentally on alto saxophone fingerings *Corresponding author: joel.gilbert@univ-lemans.fr during a project for the design of microinterval systems [9]. An decreased harmonicity by extending the bore of the Tintignac carnyx improves its ease of playing [10,11]. The harmonicity of resonances is also necessary for proper intonation when a reed instrument is expected to play in upper registers, and is often used as a target in optimisation problems [12][13][14][15][16]. This paper focuses on the assessment of the Bouasse-Benade prescription on a model of reed musical instruments. The goal is to study the influence of the inharmonicity on the playing frequency and on the minimum mouth pressure required to achieve auto-oscillation in the first register. However, it appears that this mouth pressure cannot always be determined by a study of small amplitude oscillations only. On the contrary, a complete bifurcation diagram, including all periodic branches with the blowing pressure as the continuation parameter, needs to be computed. To achieve this, a simplified model of a reed instrument is derived from a generic model that is valid for both reed and brass instruments, and constitutes a simplified version of the problem. The reed is modeled as a simple spring [17], only two acoustic resonances are taken into account [18,19], and the nonlinear coupling between the reed and the acoustic resonances through the incoming flow is reduced to a polynomial expansion (Kergomard in [20]). This may be considered as the simplest model of reed instruments that includes inharmonicity. Furthermore, this simple model helps isolate the effects of the main parameters without the added complications that arise when considering real instruments.
Inharmonicity Inh between the acoustic resonances f res1 and f res2 is defined as the deviation from harmonicity: Inh = f res2 /(2f res1 ) À 1. Therefore once the first resonance frequency f res1 and the Inharmonicity Inh are known, the second resonance frequency is fixed through the relation f res2 = 2f res1 (1 + Inh). Note that if the resonances are exactly harmonic (Inh = 0) the problem can be solved analytically and two bifurcation diagrams have been obtained (Figs. 8 and 10 of [18]).
In Section 2 of this paper the theoretical background, and particularly the equations of the elementary model of reed instruments, are briefly presented. The behaviour of the elementary model at the stability threshold of the equilibrium position, and the nature of the Hopf bifurcations, are discussed in the second part of this section. Section 3 documents the procedure used to calculate bifurcation diagrams using a continuation method, after having reformulated the two equations of the model into a set of four first-order ODE equations. The influence of the inharmonicity on the bifurcation diagrams is shown and discussed in Section 4. The section is divided into two parts: the first study assumes that the amplitude of the first resonance is larger than the second resonance (Z 1 > Z 2 ), for which preliminary results have been presented in [21], and the second study the opposite condition (Z 2 > Z 1 ) is considered. In order to link the Bouasse-Benade prescription to the ease of playing experienced by musicians, bifurcation diagrams are analysed with respect to the minimal mouth pressure necessary to achieve oscillation. Additionally, the effect of inharmonicity is also considered. The model presented and used in the present publication is labeled as elementary because a number of major simplifications are made in deriving it (see for example Hirschberg in [20,22]). The vibrating reeds or lips are modeled as a linear one-degree-of-freedom oscillator. The upstream resonances of the player's windway are neglected, as is nonlinear propagation of sound in the air column of the instrument. Wall vibrations are also ignored. Despite these simplifications, the elementary model is capable of reproducing many of the important aspects of performances by human players on realistic reed and brass instruments (see [4,5]). The model is based on a set of three equations, which have to be solved simultaneously to predict the nature of the sound radiated by the instrument. These three constituent equations of the model are presented hereafter. Besides the control parameters defining the embouchure of the players, including the reed or lips parameters and the mouth pressure P m , and the input impedance of the wind instrument, there are three variables in the set of the forthcoming three equations as a function of the time t:hðtÞ, the reed or lip-opening height,pðtÞ the pressure in the mouthpiece of the instrument, andũðtÞ the volume flow entering the instrument.
In order to describe the vibrating reeds or lips, the first of the three constituent equations of the elementary model is: In this equation, which describes the reeds or lips as a onedegree-of-freedom (1DOF) mechanical oscillator, the symbols x r , Q r , h o and l represent the angular reed resonance frequency, the quality factor of the reed resonance, the value of the reed or lip-opening height at rest, and the effective mass per unit area of the reed or lips respectively. These quantities are parameters of the model, which are either constant (in a stable note) or changing slowly in a prescribed way (in a music performance). Note that if l is positive, an increase of the pressure difference ðP m ÀpðtÞÞ will imply a closing of the reed or lips aperture. It is called the "inward-striking" model, used mainly for reed instruments. If l is negative, an increase of the pressure difference will imply an opening of the reed or lips aperture. It is called the "outward striking" model, used preferably for brass (lip reed) instruments. The second constituent equation describes the relationship between pressure and flow velocity in the reed channel: where the square root originates from the Bernoulli equation, and the positive part of the reed or lips aperturẽ h þ ¼ maxðh; 0Þ implies that the volume flow vanishes when the reed or lips are closed.
The third and last constituent equation describes the relationship between flow and pressure in the instrument mouthpiece. It is written in its frequency domain form by using the input impedance Z(x) of the wind instrument: Other than the difference of sign of l between inward-striking reed instruments model and outward striking brass instruments model, there is another difference between these two subfamilies of wind instruments. The control parameter x r of vibrating lips varies a lot, over four octaves, to get the entire tessitura of a given brass instrument. On the other hand the x r associated to reeds is more fixed (slightly varying because of the lower lip of the clarinet or saxophone player) and most of the time very large compared to the playing frequencies. This justifies a low-frequency approximation of the elementary model: x r is assumed infinite and the reed undamped. In other words, the reed is reduced to its stiffness only and the set of three equations becomes a set of two equations as follows: When the mouth pressure is too high, the reed can be blocked against the lay of the mouthpiece. Then the closure pressure defined by P M ¼ lx 2 r Á h o is the minimal mouth pressure for which the reed remains closed in the static regime (h becomes equal to 0). By using this closure pressure, a dimensionless mouth pressure c can be defined: c = P m /P M . It is convenient to define another dimensionless parameter, a dimensionless reed height at rest: where Z c ¼ qc S is the characteristic impedance for plane wave inside the resonator of input cross section S, q is the air density and c is the sound velocity.
In the following, the nonlinear equation of the model is approximated by its third-order Taylor series around the equilibrium position defined byp eq ¼ 0, . The approximated nonlinear equation can be written in the following dimensionless form (see for example Kergomard in [20]): p . The value of the dimensionless reed height at rest f is chosen to be equal to 0:1. It is this elementary low-frequency model for reed instruments which is studied in the present paper. If a non-beating reed is assumed which is typically obtained for a dimensionless mouth pressure c lower than 0.5, the third order approximation of the flow rates is appropriate. The elementary model based on the set of two equations has to be solved to predict the nature of the sound radiated by the instrument. Low amplitude solutions for a few specific cases are reviewed in the following subsection.

Small amplitude behaviour
The equilibrium position is the trivial permanent (steady) regime corresponding to silence. Sound can happen if the equilibrium position becomes unstable. For a lossless cylindrical air column, it becomes unstable for a specific value of c which is c thr = 1/3. If losses are taken into account, then c thr is a bit higher (see [5]). If the losses are very important, the threshold value c thr can reach 1 and the reed channel is closed at equilibrium. In this case the equilibrium remains stable for any value of c. Hence no sound can be produced. For an extensive analysis of stability of the equilibrium position with an experimental comparison for cylindrical air columns, see [23] and [24].
The step beyond the above linear stability analysis is the study of the small oscillations around the threshold. It has been done first by [25] and then extended by analysing the nature of the bifurcation at the threshold by [26] which can be direct or inverse Hopf bifurcation. The results are displayed in Figure 1 as a 2D map where the x-axis is C the third coefficient of the Taylor expansion Equation (6), and the y-axis is 1/Z 2 À 1/Z 1 , the difference between the admittance amplitude between the two first resonances (assumed to be harmonic, the ratio between their frequencies, being equal to 2).
In our specific case, the coefficient C is negative. Then, for a specific negative C value, following an imaginary vertical line coming from an infinite positive value of 1/Z 2 À 1/Z 1 (second resonance peak absent like for the cylindrical tube) the bifurcation is direct. It becomes inverse in a particular point for a particular positive value of 1/Z 2 À 1/Z 1 not far from zero: 1/Z 2 À 1/Z 1 = À2B 2 / (3C). And when 1/Z 2 À 1/Z 1 becomes negative, and whatever how 1/Z 2 < 1/Z 1 (it means whatever Z 2 > Z 1 ) is, the bifurcation becomes and stays direct. Properties of small amplitude oscillations of the single-reed woodwind instruments near the oscillation threshold have been investigated more recently by using analytical formulae with explicit dependence on the physical parameters of the instrument and the instrumentist allowing to determine the bifurcation point, the nature of the bifurcation, the amplitude of the first harmonics and the oscillation frequency [27]. Apart from a few very simplified cases, such as a clarinet-like model with a lossless cylindrical tube ( [28], Kergomard in [20,5,29]), or by taking into account losses independent of frequency, sometimes called Raman model [30], the equations are not tractable analytically, and the bifurcation diagrams can not be easilly obtained.
The simplest non-trivial resonator that can be studied, is a resonator having two quasi-harmonic resonance frequencies f res1 and f res2 . This kind of resonator can be obtained in practice in the midle and high ranges of the first register of saxophone (see for example Figs. 17 and 12 of [9]). Bifurcation diagrams have been analytically calculated in [18] in the restrictive case of perfect harmonicity between the two resonances. In the following sections, this kind of resonator but with a non-zero inharmonicty Inh is analysed.

Typical bifurcation diagram obtained by continuation method
To overcome the difficulties of the analytical analysis of small amplitude oscillations near thresholds, and to get results for any inharmonicity value arbitrarily far from the oscillation threshold, simulation techniques in time domain are often used. An alternative method is possible. A nice way to have an overview of the dynamics over small and large amplitudes is to use the bifurcation diagram representation. Very few of them can be obtained analytically (see the previous subsection). It is possible to obtain bifurcation diagrams numerically for a large range of situations by using continuation methods, such as in the AUTO software [31] or MANLAB software [32] for example. In order to use AUTO technique in the following section, the elementary model has to be mathematically reformulated in a set of first-order ODE equations.
The principle of continuation is to seek solution branches of a nonlinear algebraic system rather than solution points. A solution branch is an 1D-curve in a space whose axes are an unknown to the problem and a parameter of interest called a bifurcation parameter. In the following the dimensionless mouth pressure c is chosen as the bifurcation parameter. It provides more information than a set of solution points obtained for successive values of the bifurcation parameter. Branches of static and periodic solutions are computed numerically hereafter using the software AUTO, freely available online [33].
The model analysed in this paper is a nonlinear dynamical system. In order to obtain a nonlinear algebraic system in which numerical continuation can be applied, some additional work may be required. For instance, for continuing periodic solutions of a dynamical system, a discretisation is necessary to come down to an algebraic system. Many approaches are possible among which a time-domain discretisation of the (unknown) solution over one (unknown) period. The unknowns of the resulting nonlinear algebraic system are the sampled values of the periodic solution and the period. The time discretisation implemented in AUTO is called orthogonal collocation and relies on the use of Lagrange polynomials. The stability of each solution is also assessed. Stability is a very important information for the interpretation of the bifurcation diagram since only stable solutions are observable. Stability of both equilibria and periodic solutions is found through a linearization of the system of equations around the solution considered. The solution is stable if and only if the real parts of all the eigenvalues of a matrix characteristic of the linearized system are negative. This matrix is the Jacobian matrix if the solution considered is an equilibrium, and the so-called monodromy matrix if the solution considered is periodic. Stability of a solution along a branch is an output of AUTO. For comprehensive details about continuation of static/periodic solutions using AUTO, please refer to [31].
In order to use the AUTO technique, the input impedance equation (Eq. (3)) is reformulated by a sum of individual acoustical resonance modes in the frequency domain, and then translated them in the time domain. There are two ways to manage that: sum of real modes (see for example [34]), sum of complex modes (see for example [35]). These two ways of approximating the input impedance in the frequency domain lead to two different sets of first-order equations dX dt ¼ F ðX Þ with two different X vectors. In the present paper the real mode representation of the input impedance Z is used.
The modal-fitted input impedance with N resonance modes, is written as follows: where the nth resonance is defined by three real constants, the amplitude Z n , the dimensionless quality factor Q n and the angular frequency x n . Translation of Equation (7) in the time domain and reconstruction of p(t) from real modal components p n , such that the acoustical pressure is pðtÞ ¼ P N n¼1 p n ðtÞ, results in a second order ODE for each p n : Taking into account the other equation of the elementary model, the time derivative of the volume flow nonlinear equation (Eq. (6)), the previous set of N second order  (6), and the y-axis is the difference between the admittance amplitude between the two first resonances (assumed to be harmonic) Y 2 ÀY 1 = 1/Z 2 À 1/Z 1 . The hatched region is for a direct bifurcation, and the unhatched region for an inverse bifurcation. Adapted from [26].
ODE (Eq. (8)) can be rewritten by using the following expression of du dt : Then the equations can be put into a state-space representation dX dt ¼ F ðX Þ, where F is a nonlinear vector function, and X the state vector having 2N real components defined as follows: . . . ; p N ; dp 1 dt ; . . . ; dp N dt In practice, because our paper is dedicated to a two quasi harmonic resonance instrument, the state space representation is based on the state vector of four real components X ¼ p 1 ; p 2 ; dp 1 dt ; dp 2 dt Â Ã , and the nonlinear vector function F can be written as: Before discussing extensively bifurcation diagrams for different values of inharmonicity and for different configurations of relative amplitudes between Z 1 and Z 2 of the two resonances, let us begin by showing and discussing typical elements of a bifurcation diagram. Figure 2 has been obtained by choosing Z 1 = 1.5Z 2 and two harmonic resonances (i.e. Inh = 0). The values of the modal parameters of the two resonance's air column given in Table 1 are inspired from [9] and [19]. The main plot displays the continuation results obtained with AUTO, whereas the six smaller plots above correspond to time-domain simulations of the same system between t = 0 and t = 0.5 s, for different values of c pointed by numbers. Time integration is performed with an ordinary differential equation solver, namely ode15s from the Matlab ODE Suite. The main plot displays max|p|, the maximum of the absolute value of pressure in the mouthpiece over one period with respect to the blowing pressure c. While it is not highlighted here, the horizontal line max|p| = 0 corresponds to the equilibrium solution. Below a certain critical value of c (namely c < c thr1 ), the equilibrium is stable as illustrated by the three time domain simulations calculated for c = 0.32, c = 0.36 and c = 0.40. For initial conditions chosen around the equilibrium, these oscillating solutions decay in time back to the (stable) equilibrium. It is worth noting that the decaying transient lasts all the longer as the value of c is approaching the critical value c thr1 . When c = c thr1 , the equilibrium becomes unstable and a branch of periodic solution emerges from the equilibrium. This branch is represented in green on the main plot of Figure 2: it first goes backward in terms of c and is unstable (thin line), then after a turning point (also called a fold) goes forward and is stable (thick line). This scenario is called an inverse Hopf bifurcation and the value c = c subthr the sub-critical threshold (see, e.g. [36]).
As explained above, the bifurcation point c = c thr1 is reached when the real part of one eigenvalue of the jacobian matrix crosses the imaginary axis. The imaginary part of the eigenvalue concerned gives the angular frequency of the emerging periodic solution. In the present case, it is close to x 1 . Hence the periodic solution is classified as "first register" or fundamental regime. If the angular frequency of the emerging periodic solution had been close to x 2 , the periodic solution would have been classified as "second register" or octave regime. Note that the frequency of the periodic solution along the green branch is not locked at any value but is modified according to the nonlinearity. This is exemplified and discussed in the next section. Two time domain simulations are shown with c = 0.4 and c = 0.45 and reveal that the solution is repelled from the equilibrium and converges toward a periodic solution. Note that in the case of c = 0.4 the choice of the initial condition is crucial since two stable solutions exist: the equilibrium (plot number 3 in Fig. 2) and the periodic solution (plot number 4). A thorough look at the time domain simulation would reveal that max|p| deduced from the steady-state (periodic) regime is equal to the ordinate of the green curve at the corresponding value of c.
The black curve corresponds to emerging branch of periodic solutions in the case where only one acoustic resonance is considered (Z 2 = 0). In that case, the amplitude max|p| is simply a square-root shaped function of the bifurcation parameter c in the neighbourhood of the threshold. The thick line denotes a stable periodic solution. Such a scenario is called a direct Hopf bifurcation. Just above the Hopf bifurcation point (c = c thr1 ), the direct scenario leads to stable periodic oscillations with infinitely small amplitudes. Sounds can be played with the nuance pianissimo. On the contrary, in the case of an inverse bifurcation, stable periodic oscillations found just above the Hopf bifurcation point have finite amplitude. Playing with the pianissimo nuance is no longer possible.
For pedagogical purposes, the bifurcation diagram is limited here to the neighborhood of one Hopf bifurcation point, coming from the value c = c thr1 . However, it will be shown in the next section that for other values of c, another Hopf bifurcation point is found as well as other bifurcations of the periodic branches.

Effects of the inharmonicity
Results and discussion

Large first resonance amplitude
The discussion is initiated by analysing the case corresponding to Z 1 slightly higher than Z 2 (in practice Z 1 /Z 2 = 3/2). Three bifurcation diagrams corresponding to Inh = 0, Inh = 0.02 and Inh = 0.04 are shown in Figure 3 (remember that a semi tone corresponds to 0:059).   The results shown in Figure 3 for the case Inh = 0 are qualitatively consistent with the one published in [18] (see in particular its Fig. 8). Note that the continuation method gives an additional information: the stability nature of the periodic oscillations.
In Figure 3a the bifurcation diagram shows two branches coming from the equilibrium position: 1. The first branch originates from the linear threshold c = c thr1 , associated to the first resonance f res1 , originating through an inverse Hopf bifurcation. This fundamental regime, or first register regime, is a standard Helmholtz motion according to [18]. The branch is unstable and then becomes stable at the limit point at c = c subthr (sub-critical threshold). Compared to the case of a single mode (black curve), important differences are observed, including the nature of the bifurcation. 2. The second branch originates from the linear threshold c = c thr2 and is associated to the second resonance f res2 , originating through a direct Hopf bifurcation. Note that c thr2 is above c thr1 , because Z 1 is bigger than Z 2 . This branch which would correspond to the octave regime, or second register regime, is not observable in practice, because the periodic solutions are unstable. 3. The nature of the bifurcation of the two branches originating from the linear thresholds c = c thr1 and c = c thr2 is in agreement with the publication of [26]. 4. There is a third branch which originates from the unstable octave branch, thanks to a period doubling bifurcation. This branch which would correspond to another fundamental regime (the inverted Helmholtz motion according to [18]) is unstable.
The associated lower plot shows the frequency of the periodic oscillations corresponding to the branches of the bifurcation diagram. In particular the frequency of the fundamental regime (green curve) is almost locked to the value f res1 = f res2 /2 for any value of c.
For an inharmonicity of 0.02 (Fig. 3b). The bifurcation diagram is quite close to the one with Inh = 0. However two things are pointed out. First, at the threshold c = c thr1 the Hopf bifurcation has become direct as it can be predicted theoretically [37]. Second, again there are periodic oscillations for values of the mouth pressure c under c = c thr1 until a new value c = c subthr which is a bit larger than the one of the case Inh = 0. This is due to the occurence of two folds (limit points on the solution branch) corresponding to saddle-node bifurcations. Note that the frequency of the fundamental regime (green curve) is not locked at the value f res1 anymore but is partially pulled towards the value f res2 /2, which is reasonable. If the inharmonicity was negative, the same kind of results would have been obtained, the frequency being pulled towards f res2 /2 lower than f res1 .
For an inharmonicity of 0.04 (Fig. 3c). Now the branch coming from the threshold c = c thr1 , corresponding to the fundamental regime, looks like a classical branch associated to the direct Hopf bifurcation, there is no c subthr anymore, since the folds noted in the previous case have disappeared, c = c thr1 is now the threshold of oscillation. In fact, when the inharmonicity increases, the dynamics of the system behaves more and more like the dynamics of a singleresonance system. The frequency of the fundamental regime comes from the threshold value f thr1 at the direct Hopf bifurcation point, and then is partially pulled toward the value f res2 /2. Note that in Figure 3a curve corresponding to a direct Hopf bifurcation is branched at c = c thr1 , this curve corresponds to an air column having only one resonance at the frequency f res1 .
Under certain circumstances, for instance when the inharmonicity is high enough, a branch of quasi-periodic solutions may emerge from a Neimark-Sacker bifurcation (often refered as a Hopf bifurcation for a periodic regime). Above this bifurcation point, the periodic branch still exists but it becomes unstable. Such bifurcation has not been encountered in this work, but it has been observed experimentally with a modified saxophone played in the medium range of its tessitura [9], simulated by [38], extensively studied in [19], and it has been studied with continuation on a toy model of saxophone in Section 3 of [39].
The above analysis illustrates significant things because of the inverse Hopf bifurcation (cases Inh = 0 and Inh = 0.02): 1. On the one hand, there may be a minimum value c = c subthr lower than c thr1 above which there are stable periodic oscillations. This particular value c subthr can be thought of as a quantitative characterisation of the ease of playing. In Figure 3 it is shown that the lowest value of c subthr is obtained when the two resonances are perfectly harmonic (Inh = 0). If it is assumed that a lower c subthr corresponds to an instrument easier to play, then it suggests the reed instrument considered is the easiest to play when Inh = 0. In a way that is a theoretical illustration of the Bouasse-Benade prescription. The threshold of oscillation, equal to c subthr for low inharmonicities, and equal to c thr1 for higher inharmonicities, is displayed in Figure 4. The minimum of the threshold of oscillation correspond to Inh = 0. 2. On the other hand, the stable periodic oscillations which appear for c slightly above c subthr can have fundamental frequencies significantly different from f thr1 = f res1 because of the effect of the second resonance which controls partially the intonation of the fundamental regime. This study highlights the intrinsic limitation of the linear stability analysis: it should be only considered to assess the stability of the equilibrium. Conclusions concerning the existence of periodic solutions can only be provided through a nonlinear analysis, analitycally in specific cases or with tools like AUTO otherwise.
In addition an animation showing the evolution of the bifurcation diagram as a function of the inharmonicity increasing from Inh = À0.05 to Inh = +0.05 is available from the link of footnote 1 . Most of the illustrations displayed in the figures are corresponding to a positive inharmonicity Inh, but the animation and the Figure 4 illustrate the fact that the behaviour is qualitavely the same for negative values of Inh.
In order to illustrate the bifurcation diagram (Fig. 3), it is interesting to do simulations by solving the equation dX dt ¼ F ðX Þ in the time domain (sounds available from the links of footnotes 2, 3 ). Figures 5 and 6 show a signal corresponding to an inharmonicity Inh = 0.040 and Z 1 = 1.5Z 2 : 1. In Figure 5 the control parameter increases linearly from c = 0.43 to c = 0.50 (crescendo). Because the branch is coming from a direct Hopf bifurcation in the bifurcation diagram, the amplitude of the oscillation is a smoothly increasing mathematical function with respect to c. Therefore, in the time domain simulation, the amplitude of the signal (fundamental regime) is increasing smoothly, as it is with a resonator having only one resonance f res1 . 2. In Figure 6 the control parameter decreases slowly from c = 0.55 to c = 0.50 (decrescendo). Because of the chosen initial conditions, the periodic regime obtained is corresponding to the upper octave, but when c reaches the value 0.53, the branch coming from c = c thr2 becomes unstable, and then the periodic solution jumps on the first branch one octave below, the stable branch coming from c = c thr1 (fundamental regime).

Large second resonance amplitude
The discussion continues by analysing the case corresponding to Z 1 slightly lower than Z 2 (in practice Z 2 /Z 1 = 3/2). Three bifurcation diagrams corresponding to Inh = 0, Inh = 0.015 and Inh = 0.03 are shown Figure 7.
The results shown in Figure 7 for the case Inh = 0 are qualitatively consistent with the one published in [18] (see in particular its Fig. 10).
In Figure 7a the bifurcation diagram shows two branches coming from the equilibrium position: 1. On the left-hand side, the first branch originates from the linear threshold c = c thr2 , associated to the second resonance f res2 , originating through a direct bifurcation. This octave regime is stable until a period doubling bifurcation point, and then becomes unstable. At the bifurcation point, there is an emerging branch corresponding to a fundamental regime. It is a standard Helmholtz motion according to [18]. This fundamental regime is unstable until a turning point (a fold) corresponding to a minimum value of c = c subthr where the periodic oscillations become stable. Note that the threshold of oscillation of the fundamental regime c = c subthr is significantly lower than the value c thr1 predicted by the linear stability analysis.    Fig. 3c). The dimensionless mouth pressure is printed in black, and is decreasing linearly from c = 0.55 (constant before t = 2 s) to c = 0.50 (constant after t = 9 s).
2. The second branch originates from the linear threshold c = c thr1 , associated to the first resonance f res1 , originating through a direct bifurcation. Note that c thr1 is bigger than c thr2 , because Z 1 is lower than Z 2 . This branch which would correspond to a second fundamental regime is not observable in practice, because the periodic solutions are unstable. This branch would correspond to the inverted Helmholtz motion according to [18].
The associated lower curve shows the frequency of the periodic oscillations corresponding to the branches of the bifurcation diagram. In particular the frequency of the stable fundamental regime (blue curve) is close to the value f res1 = f res2 /2 for any value of c.
For an inharmonicity of 0.015 (Fig. 7b). The bifurcation diagram is qualitatively quite close to the one with Inh = 0. Two things are now pointed out. Once again there are periodic oscillations for dimensionless mouth pressure c values below c = c thr2 < c thr1 until a new value c = c subthr which is a bit bigger than the one in the case of Inh = 0. Note that the frequency of the fundamental regime (blue curve) is surprisingly close to the value f res1 , the fundamental frequency is not much pulled towards the value f res2 /2. Note that the range of c where there are two stable periodic regimes, octave and the standard Helmholtz motion fundamental regime, is larger: from c = c thr2 to the value of c where the period doubling bifurcation point occurs.
For an inharmonicity of 0.03 (Fig. 7c). Again the bifurcation diagram is qualitatively quite close to the ones with Inh = 0 and Inh = 0.015. The minimum pressure of fundamental periodic oscillations c subthr (on the blue branch keeps increasing with inharmonicity, and becomes higher than c = c thr2 . The above discussion illustrates significant things: 1. There may be a minimum value c = c subthr lower than c thr2 < c thr1 where there are stable periodic oscillations. This particular value c subthr can be chosen as a kind of quantitative characterisation of the ease of playing. In Figure 8 it is shown again (as in Fig. 4) that the lowest value of c subthr is obtained when the two resonances are perfectly harmonic (Inh = 0). If it is assumed that a lower c subthr corresponds to an instrument easier to play, then it suggests that the reed instrument considered to be the easiest to play when Inh = 0. In a way, even if Z 1 < Z 1 , again that is a theoretical illustration of the Bouasse-Benade  prescription. The threshold of oscillation is displayed in Figure 8: the minimum of the threshold of oscillation is corresponding to Inh = 0. 2. The stable periodic oscillations which appear for c slightly bigger than c subthr have fundamental frequencies quite close to f res1 . 3. It is worth emphasising that, whatever the inharmonicity, the fundamental regime does never come from the first threshold c = c thr1 , but comes through a period-doubling bifurcation point attached to the octave branch. A naive analysis of the time domain simulations (at least with Inh = 0) would probably suggest that the fundamental regime emerges from the equilibrium trough an inverse Hopf bifurcation, but this is not correct. It is also worth noting that a linear stability analysis (LSA) of the equilibrium is useless here to give some hints about the oscillation behaviour of the model. 4. Note as well that sometimes, there are several stable regimes (equilibrium position and periodic regime, or two periodic regimes) for a given value of c. For such cases, the stable regime reached is the consequence of the initial conditions.
In addition, an animation showing the evolution of the bifurcation diagram as a function of the inharmonicity increasing from Inh = À0.05 to Inh = +0.05 is available from the link of footnote 4 . Most of the illustrations displayed in the figures correspond to a positive inharmonicity Inh, but the animation and the Figure 8 illustrate the fact that the behaviour is qualitavely the same for negative values of Inh. Unlike Figure 4, it can be noted that the plot is slightly asymetric with respect to the vertical axis Inh = 0.
In order to illustrate the bifurcation diagrams (Fig. 7), it is interesting to do simulations by solving the equation dX dt ¼ F ðX Þ in the time domain (sound available from the link of footnote 5 ). Figure 9 shows a signal corresponding to an inharmonicity Inh = 0.015 and Z 1 = Z 2 /1.5. The control parameter increases slowly from c = 0.38 (just below the period doubling bifurcation) to c = 0.45 (crescendo). Because of the chosen initial conditions, the periodic regime obtained corresponds to the octave, but when c reaches the value 0.39, the branch coming from c = c th2 becomes unstable, and then the periodic solution jumps to the only stable branch one octave below (branch coming from the period-doubling bifurcation).

Conclusion
Bifurcation diagrams of a basic reed instrument modeled by two quasi-harmonic resonances have been computed by using a continuation method (AUTO software), where the mouth pressure is the control parameter. Some of the mouth pressure thresholds results are interpreted in terms of the ease of playing of the reed instrument. When there is an inverse Hopf bifurcation (perfect harmonicity) or of a double fold after a direct Hopf bifurcation (moderate inharmonicity), there may be a minimum value c = c subthr lower than c thr1 for which periodic stable oscillations can be observed. This value c subthr may be considered as a quantitative characterisation of the ease of playing. It has been shown that the lowest value of c subthr is obtained when the two resonances are harmonic, harmonicity equal to 2. This is a theoretical illustration of the Bouasse-Benade prescription [2,8]). Even if a few AUTO simulations using other parameter's values than the one used in the present study have been done, a large set of other tests should be done with many other parameter's values to verify that the conlusions of the present paper are robust.
An interesting direction for future work could include experimental validation, particularly using the modified saxophone used in [9]. There, the saxophone was modified by the addition of two closed side tubes on the neck. Movable pistons are used to change the volume of the side tubes, which results in a shift of the resonant frequencies.
As explained in this publication, it is possible to choose particular positions and volumes of the closed tubes to ensure a control of the inharmonicity. This strategy is used with fingerings corresponding to the middle and high ranges of the first regime of the saxophone, where its input impedance consists essentially of two resonances.
The results provided in the current manuscript depend on a physical model of reed instruments based on three strong approximations: the reed dynamics is ignored, only two acoustic resonances are taken into account, and the nonlinear equation describing the incoming volume flow is approximated by its third order Taylor series expansion. Therefore, the conclusions above cannot be directly extended to real instruments until further research is carried out on more complex models. At that point, many other interesting topics could be explored, such as sound  Fig. 7b). The dimensionless mouth pressure is printed in black, and is increasing linearly from c = 0.38 (constant before t = 2 s) to c = 0.45 (constant after t = 9 s). production of low-pitched notes by conical reed instruments such as saxophones, oboes or bassoons. Replacing the inward-striking reed model by an outward-striking lip model is also planned in order to study nonlinear dynamics of brass instruments (preliminary results in [40]). More precisely, it is expected that numerical continuation could clarify their pedal note regime, recently simulated using a time-domain finite differences method in [41].