© J. Gilbert et al., Published by EDP Sciences, 2020

1 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–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 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–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₁ > Z₂), for which preliminary results have been presented in [21], and the second study the opposite condition (Z₂ > Z₁) 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.

2 Theoretical background

2.1 Elementary acoustical model

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: $\tilde{h}\left(t\right)$, the reed or lip-opening height, $\tilde{p}\left(t\right)$ the pressure in the mouthpiece of the instrument, and $\tilde{u}\left(t\right)$ 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:$$\frac{{\mathrm{d}}^2\tilde{h}\left(t\right)}{\mathrm{d}{t}^2}+\frac{{\omega }_r}{Q_r}\frac{\mathrm{d}\tilde{h}\left(t\right)}{\mathrm{d}t}+{\omega }_r^2\left(\tilde{h}\left(t\right)-h_o\right)=-\frac{P_{\mathrm{m}}-\tilde{p}\left(t\right)}{\mu }.$$(1)

In this equation, which describes the reeds or lips as a one-degree-of-freedom (1DOF) mechanical oscillator, the symbols ω_r, Q_r, h_o and μ 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 μ is positive, an increase of the pressure difference $(P_{\mathrm{m}}-\tilde{p}(t\left)\right)$ will imply a closing of the reed or lips aperture. It is called the “inward-striking” model, used mainly for reed instruments. If μ 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:$$\tilde{u}\left(t\right)=w{\tilde{h}}^{+}\left(t\right)\sqrt{\frac{2}{\rho }\left|P_m-\tilde{p}\left(t\right)\right|}\mathrm{sign}\left(P_m-\tilde{p}\left(t\right)\right),$$(2)where the square root originates from the Bernoulli equation, and the positive part of the reed or lips aperture ${\tilde{h}}^{+}=\max (\tilde{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(ω) of the wind instrument:$$\tilde{P}\left(\omega \right)=Z\left(\omega \right)\tilde{U}\left(\omega \right).$$(3)

Other than the difference of sign of μ 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 ω_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 ω_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: ω_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:$$\{\begin{array}{l}\tilde{u}\left(t\right)=w\left[h_{\mathrm{o}}-\frac{P_{\mathrm{m}}-\tilde{p}\left(t\right)}{\mu {\omega }_r^2}\right]\sqrt{\frac{2}{\rho }|P_{\mathrm{m}}-\tilde{p}\left(t\right)|}\mathrm{sgn}\left(|P_{\mathrm{m}}-\tilde{p}\left(t\right)|\right)\\ \tilde{P}\left(\omega \right)=Z\left(\omega \right)\tilde{U}\left(\omega \right).\end{array}$$(4)

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=\mu {\omega }_r^2·h_o$ is the minimal mouth pressure for which the reed remains closed in the static regime ($\tilde{h}$ becomes equal to 0). By using this closure pressure, a dimensionless mouth pressure γ can be defined: γ = P_m/P_M. It is convenient to define another dimensionless parameter, a dimensionless reed height at rest:$$\zeta =Z_{c}w{h}_o\sqrt{\frac{2}{\rho P_M}},$$(5)where $Z_c=\frac{\rho c}{S}$ is the characteristic impedance for plane wave inside the resonator of input cross section S, ρ 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 by ${\tilde{p}}_{\mathrm{eq}}=0$, ${\tilde{h}}_{\mathrm{eq}}=h_o-\frac{P_{\mathrm{m}}}{\mu {\omega }_r^2}$, and ${\tilde{u}}_{\mathrm{eq}}=w{\tilde{h}}_{\mathrm{eq}}\sqrt{\frac{2}{\rho }{P}_{\mathrm{m}}}$.

The approximated nonlinear equation can be written in the following dimensionless form (see for example Kergomard in [20]):$$u\left(t\right)=u_{\mathrm{eq}}+\mathrm{Ap}\left(t\right)+\mathrm{Bp}(t{)}^2+\mathrm{Cp}(t{)}^3,$$(6)with $u_{\mathrm{eq}}=\zeta (1-\gamma )\sqrt{\gamma }$, $A=\zeta \frac{3\gamma -1}{2\sqrt{\gamma }}$, $B=-\zeta \frac{3\gamma +1}{8\sqrt{{\gamma }^3}}$ and $C=-\zeta \frac{\gamma +1}{16\sqrt{{\gamma }^5}}$. The value of the dimensionless reed height at rest ζ 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 γ 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.

2.2 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 γ which is γ_thr = 1/3. If losses are taken into account, then γ_thr is a bit higher (see [5]). If the losses are very important, the threshold value γ_thr can reach 1 and the reed channel is closed at equilibrium. In this case the equilibrium remains stable for any value of γ. 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₂ − 1/Z₁, the difference between the admittance amplitude between the two first resonances (assumed to be harmonic, the ratio between their frequencies, being equal to 2).

Figure 1 Diagram showing the regions where the bifurcation is direct as well as those regions where it is inverse. The x-axis shows the values of the third coefficient of the Taylor expansion Equation (6), and the y-axis is the difference between the admittance amplitude between the two first resonances (assumed to be harmonic) Y2 −Y1 = 1/Z2 − 1/Z1. The hatched region is for a direct bifurcation, and the unhatched region for an inverse bifurcation. Adapted from [26].

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₂ − 1/Z₁ (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₂ − 1/Z₁ not far from zero: 1/Z₂ − 1/Z₁ = −2B²/(3C). And when 1/Z₂ − 1/Z₁ becomes negative, and whatever how 1/Z₂ < 1/Z₁ (it means whatever Z₂ > Z₁) 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.

3 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 γ 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 $\frac{\mathrm{d}X}{\mathrm{d}t}=F\left(X\right)$ 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:$$Z\left(\omega \right)=\sum _{n=1}^N\mathrm{}{Z}_n\frac{\mathrm{j\omega }{\omega }_n/Q_n}{{\omega }_n^2+\mathrm{j\omega }{\omega }_n/Q_n-{\omega }^2},$$(7)where the nth resonance is defined by three real constants, the amplitude Z_n, the dimensionless quality factor Q_n and the angular frequency ω_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\left(t\right)={\sum }_{n=1}^N\mathrm{}{p}_n\left(t\right)$, results in a second order ODE for each p_n:$$\frac{{\mathrm{d}}^2{p}_n}{\mathrm{d}{t}^2}+{\omega }_n/Q_n\frac{\mathrm{d}{p}_n}{\mathrm{d}t}+{\omega }_n^2{p}_n\left(t\right)=Z_n{\omega }_n/Q_n\frac{\mathrm{d}u}{\mathrm{d}t}.$$(8)

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 ODE (Eq. (8)) can be rewritten by using the following expression of $\frac{\mathrm{d}u}{\mathrm{d}t}$:$$\frac{\mathrm{d}u}{\mathrm{d}t}=A\left(\sum _{n=1}^N\mathrm{}\frac{\mathrm{d}{p}_n}{\mathrm{d}t}\right)+2B\left(\sum _{n=1}^N\mathrm{}\frac{\mathrm{d}{p}_n}{\mathrm{d}t}\right)\left(\sum _{n=1}^N\mathrm{}{p}_n\left(t\right)\right)+3C\left(\sum _{n=1}^N\mathrm{}\frac{\mathrm{d}{p}_n}{\mathrm{d}t}\right){\left(\sum _{n=1}^N\mathrm{}{p}_n\left(t\right)\right)}^2.$$(9)

Then the equations can be put into a state-space representation $\frac{\mathrm{d}X}{\mathrm{d}t}=F\left(X\right)$, where F is a nonlinear vector function, and X the state vector having 2N real components defined as follows:$$X={\left[p_1;\dots ;p_N;\frac{\mathrm{d}{p}_1}{\mathrm{d}t};\dots ;\frac{\mathrm{d}{p}_N}{\mathrm{d}t}\right]}^{\mathrm{\prime }}.$$(10)

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=\left[p_1;p_2;\frac{\mathrm{d}{p}_1}{\mathrm{d}t};\frac{\mathrm{d}{p}_2}{\mathrm{d}t}\right]$, and the nonlinear vector function F can be written as:$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{l}X\left(1\right)\\ X\left(2\right)\\ X\left(3\right)\\ X\left(4\right)\\ \end{array}\right)=\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{l}{p}_1\\ p_2\\ \frac{\mathrm{d}{p}_1}{\mathrm{d}t}\\ \frac{\mathrm{d}{p}_2}{\mathrm{d}t}\\ \end{array}\right)$$$$=\left(\begin{array}{c}X\left(3\right)\\ X\left(4\right)\\ \begin{array}{c}-\frac{{\omega }_1}{Q_1}X\left(3\right)-{\omega }_1^{2}X\left(1\right)+Z_1\frac{{\omega }_1}{Q_1}\frac{\mathrm{d}u}{\mathrm{d}t}\\ -\frac{{\omega }_1}{Q_2}X\left(4\right)-{\omega }_2^{2}X\left(2\right)+Z_2\frac{{\omega }_2}{Q_2}\frac{\mathrm{d}u}{\mathrm{d}t}\end{array}\end{array}\right).$$(11)

Before discussing extensively bifurcation diagrams for different values of inharmonicity and for different configurations of relative amplitudes between Z₁ and Z₂ 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.5Z₂ 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 γ pointed by numbers. Time integration is performed with an ordinary differential equation solver, namely ode15s from the Matlab ODE Suite.

Figure 2 Bifurcation diagram and time domain simulations of the two harmonic resonance air column (parameter’s values in Tab. 1 with Z2 = Z1/1.5 or with Z2 = 0, and Inh = 0) with respect to the control parameter γ. Upper plots: six time domain simulations of the dimensionless pressure p = p1 + p2 calculated between t = 0 and t = 0.5 s for γ = 0.32, γ = 0.36, γ = 0.40 (two simulations with different initial conditions), γ = 0.45 and γ = 0.53. The dimensionless pressure of the plots numbered from 1 to 3 (respectively 4 to 6) is displayed between −0.3 and +0.3 (respectively −1.2 and +1.2). Lower plot: Maximum of the absolute value of the periodic solution p over one period with respect to γ. The branch in green (respectively in black) corresponds to the case Z2 = Z1/1.5 (respectively Z2 = 0), and illustrates an inverse (respectively direct) Hopf bifurcation scenario. Stable (respectively unstable) solutions are plotted with thick (respectively thin) lines. For each scenario, the Hopf bifurcation point where the equilibrium becomes unstable, is noted γthr. In the case of an inverse bifurcation the subcrital threshold γsubthr is highlighted with a vertical dashed line.

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 γ. While it is not highlighted here, the horizontal line max|p| = 0 corresponds to the equilibrium solution. Below a certain critical value of γ (namely γ < γ_thr1), the equilibrium is stable as illustrated by the three time domain simulations calculated for γ = 0.32, γ = 0.36 and γ = 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 γ is approaching the critical value γ_thr1. When γ = γ_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 γ 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 γ = γ_subthr the sub-critical threshold (see, e.g. [36]).

As explained above, the bifurcation point γ = γ_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 ω₁. 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 ω₂, 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 γ = 0.4 and γ = 0.45 and reveal that the solution is repelled from the equilibrium and converges toward a periodic solution. Note that in the case of γ = 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 γ.

The black curve corresponds to emerging branch of periodic solutions in the case where only one acoustic resonance is considered (Z₂ = 0). In that case, the amplitude max|p| is simply a square-root shaped function of the bifurcation parameter γ 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 (γ = γ_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 γ = γ_thr1. However, it will be shown in the next section that for other values of γ, another Hopf bifurcation point is found as well as other bifurcations of the periodic branches.

4 Effects of the inharmonicity. Results and discussion

4.1 Large first resonance amplitude

The discussion is initiated by analysing the case corresponding to Z₁ slightly higher than Z₂ (in practice Z₁/Z₂ = 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).

Figure 3 Three bifurcation diagrams of the two quasi-harmonic resonance air column (with Z1 > Z2, in practice Z1/Z2 = 3/2; parameters values in Tab. 1) with respect to the control parameter γ. From top to bottom: Inh = 0 (a), Inh = 0.02 (b) and Inh = 0.04 (c). Each case is described with two plots. Upper plot: Maximum of the absolute value of the periodic solution with respect to γ. The branch in green (respectively in red, and in blue) corresponds to the fundamental regime, the standard Helmholtz motion (respectively the octave regime, and the inverted Helmholtz motion fundamental regime). Note that a black curve corresponding to a direct Hopf bifurcation is branched at γ = γthr1. This fundamental regime corresponds to an air column having only one resonance at the frequency fres1. Lower plot: Frequency with respect to γ. The frequency branch in green (respectively in red) corresponds to the fundamental frequency of the fundamental regime (respectively the octave regime, frequency divided by 2). The reference dashed horizontal lines are the reference frequencies: fres1 and fres2/2.

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 γ = γ_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 γ = γ_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 γ = γ_thr2 and is associated to the second resonance f_res2, originating through a direct Hopf bifurcation. Note that γ_thr2 is above γ_thr1, because Z₁ is bigger than Z₂. 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 γ = γ_thr1 and γ = γ_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 γ.

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 γ = γ_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 γ under γ = γ_thr1 until a new value γ = γ_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 γ = γ_thr1, corresponding to the fundamental regime, looks like a classical branch associated to the direct Hopf bifurcation, there is no γ_subthr anymore, since the folds noted in the previous case have disappeared, γ = γ_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 single-resonance 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 γ = γ_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 γ = γ_subthr lower than γ_thr1 above which there are stable periodic oscillations. This particular value γ_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 γ_subthr is obtained when the two resonances are perfectly harmonic (Inh = 0). If it is assumed that a lower γ_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 γ_subthr for low inharmonicities, and equal to γ_thr1 for higher inharmonicities, is displayed in Figure 4. The minimum of the threshold of oscillation correspond to Inh = 0.

   Figure 4 Minimum value of the mouth pressure γ (green line) corresponding to a stable periodic solution (fundamental regime) with respect to the inharmonicity Inh (case Z1 > Z2; in practice Z1/Z2 = 3/2). Linear threshold γthr1 (blue dashed line).

2. On the other hand, the stable periodic oscillations which appear for γ slightly above γ_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¹. 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 $\frac{\mathrm{d}X}{\mathrm{d}t}=F\left(X\right)$ in the time domain (sounds available from the links of footnotes², ³). Figures 5 and 6 show a signal corresponding to an inharmonicity Inh = 0.040 and Z₁ = 1.5Z₂:

1. In Figure 5 the control parameter increases linearly from γ = 0.43 to γ = 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 γ. 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 γ = 0.55 to γ = 0.50 (decrescendo). Because of the chosen initial conditions, the periodic regime obtained is corresponding to the upper octave, but when γ reaches the value 0.53, the branch coming from γ = γ_thr2 becomes unstable, and then the periodic solution jumps on the first branch one octave below, the stable branch coming from γ = γ_thr1 (fundamental regime).

Figure 5 Signal (dimensionless acoustical pressure) obtained by simulation in time domain with an inharmonicity Inh = 0.04 and Z1 = 1.5Z2 (like in Fig. 3c). The dimensionless mouth pressure is printed in black, and is increasing linearly from γ = 0.43 (constant before t = 2 s) to γ = 0.50 (constant after t = 9 s).

Figure 6 Signal (dimensionless acoustical pressure) obtained by simulation in time domain with an inharmonicity Inh = 0.04 and Z1 = 1.5Z2 (like in Fig. 3c). The dimensionless mouth pressure is printed in black, and is decreasing linearly from γ = 0.55 (constant before t = 2 s) to γ = 0.50 (constant after t = 9 s).

4.2 Large second resonance amplitude

The discussion continues by analysing the case corresponding to Z₁ slightly lower than Z₂ (in practice Z₂/Z₁ = 3/2). Three bifurcation diagrams corresponding to Inh = 0, Inh = 0.015 and Inh = 0.03 are shown Figure 7.

Figure 7 Three bifurcation diagrams of the two quasi-harmonic resonance air column (with Z1 < Z2, in practice Z2/Z1 = 3/2) with respect to the control parameter γ. From top to bottom: Inh = 0 (a), Inh = 0.015 (b) and Inh = 0.030 (c). Each case is described with two plots. Upper plot (a): Maximum of the absolute value of the periodic solution with respect to γ. The branch in red (respectively in blue, and in green) corresponds to the octave regime (respectively the fundamental regime associated to the standard Helmholtz motion, and the inverted Helmholtz motion fundamental regime). Lower plot (c): Frequency with respect to γ. The frequency branch in red (respectively in blue) corresponds to the fundamental frequency of the octave regime (respectively the fundamental regime). The frequency of the red branch has been divided by 2 for sake of clarity. The reference dashed horizontal lines are the reference frequencies: fres1 and fres2/2.

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 γ = γ_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 γ = γ_subthr where the periodic oscillations become stable. Note that the threshold of oscillation of the fundamental regime γ = γ_subthr is significantly lower than the value γ_thr1 predicted by the linear stability analysis.

2. The second branch originates from the linear threshold γ = γ_thr1, associated to the first resonance f_res1, originating through a direct bifurcation. Note that γ_thr1 is bigger than γ_thr2, because Z₁ is lower than Z₂. 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 γ.

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 γ values below γ = γ_thr2 < γ_thr1 until a new value γ = γ_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 γ where there are two stable periodic regimes, octave and the standard Helmholtz motion fundamental regime, is larger: from γ = γ_thr2 to the value of γ 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 γ_subthr (on the blue branch keeps increasing with inharmonicity, and becomes higher than γ = γ_thr2.

The above discussion illustrates significant things:

1. There may be a minimum value γ = γ_subthr lower than γ_thr2 < γ_thr1 where there are stable periodic oscillations. This particular value γ_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 γ_subthr is obtained when the two resonances are perfectly harmonic (Inh = 0). If it is assumed that a lower γ_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₁ < Z₁, 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.

   Figure 8 Minimum value of the mouth pressure γ (blue line) corresponding to a stable periodic solution (fundamental regime) with respect to the inharmonicity Inh (case Z1 < Z2; in practice Z2/Z1 = 3/2). Linear threshold γthr1 (green dashed line).

2. The stable periodic oscillations which appear for γ slightly bigger than γ_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 γ = γ_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 γ. 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⁴. 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 $\frac{\mathrm{d}X}{\mathrm{d}t}=F\left(X\right)$ in the time domain (sound available from the link of footnote⁵). Figure 9 shows a signal corresponding to an inharmonicity Inh = 0.015 and Z₁ = Z₂/1.5. The control parameter increases slowly from γ = 0.38 (just below the period doubling bifurcation) to γ = 0.45 (crescendo). Because of the chosen initial conditions, the periodic regime obtained corresponds to the octave, but when γ reaches the value 0.39, the branch coming from γ = γ_th2 becomes unstable, and then the periodic solution jumps to the only stable branch one octave below (branch coming from the period-doubling bifurcation).

Figure 9 Signal (dimensionless acoustical pressure) obtained by simulation in time domain with an inharmonicity Inh = 0.015 and Z1 = Z2/1.5 (like in Fig. 7b). The dimensionless mouth pressure is printed in black, and is increasing linearly from γ = 0.38 (constant before t = 2 s) to γ = 0.45 (constant after t = 9 s).

5 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 γ = γ_subthr lower than γ_thr1 for which periodic stable oscillations can be observed. This value γ_subthr may be considered as a quantitative characterisation of the ease of playing. It has been shown that the lowest value of γ_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 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].

Conflict of interest

The authors declare no conflict of interest.

Acknowledgments

Authors acknowledge their colleagues Murray Campbell, Jean-Pierre Dalmont and Erik Petersen for fruitfull discussions. And authors wish to thank the reviewers, whose conscientious work has improved the clarity of this article and raised some interesting questions.

References

All Tables

All Figures

Figure 1 Diagram showing the regions where the bifurcation is direct as well as those regions where it is inverse. The x-axis shows the values of the third coefficient of the Taylor expansion Equation (6), and the y-axis is the difference between the admittance amplitude between the two first resonances (assumed to be harmonic) Y2 −Y1 = 1/Z2 − 1/Z1. The hatched region is for a direct bifurcation, and the unhatched region for an inverse bifurcation. Adapted from [26].

In the text

Figure 2 Bifurcation diagram and time domain simulations of the two harmonic resonance air column (parameter’s values in Tab. 1 with Z2 = Z1/1.5 or with Z2 = 0, and Inh = 0) with respect to the control parameter γ. Upper plots: six time domain simulations of the dimensionless pressure p = p1 + p2 calculated between t = 0 and t = 0.5 s for γ = 0.32, γ = 0.36, γ = 0.40 (two simulations with different initial conditions), γ = 0.45 and γ = 0.53. The dimensionless pressure of the plots numbered from 1 to 3 (respectively 4 to 6) is displayed between −0.3 and +0.3 (respectively −1.2 and +1.2). Lower plot: Maximum of the absolute value of the periodic solution p over one period with respect to γ. The branch in green (respectively in black) corresponds to the case Z2 = Z1/1.5 (respectively Z2 = 0), and illustrates an inverse (respectively direct) Hopf bifurcation scenario. Stable (respectively unstable) solutions are plotted with thick (respectively thin) lines. For each scenario, the Hopf bifurcation point where the equilibrium becomes unstable, is noted γthr. In the case of an inverse bifurcation the subcrital threshold γsubthr is highlighted with a vertical dashed line.

In the text

Figure 3 Three bifurcation diagrams of the two quasi-harmonic resonance air column (with Z1 > Z2, in practice Z1/Z2 = 3/2; parameters values in Tab. 1) with respect to the control parameter γ. From top to bottom: Inh = 0 (a), Inh = 0.02 (b) and Inh = 0.04 (c). Each case is described with two plots. Upper plot: Maximum of the absolute value of the periodic solution with respect to γ. The branch in green (respectively in red, and in blue) corresponds to the fundamental regime, the standard Helmholtz motion (respectively the octave regime, and the inverted Helmholtz motion fundamental regime). Note that a black curve corresponding to a direct Hopf bifurcation is branched at γ = γthr1. This fundamental regime corresponds to an air column having only one resonance at the frequency fres1. Lower plot: Frequency with respect to γ. The frequency branch in green (respectively in red) corresponds to the fundamental frequency of the fundamental regime (respectively the octave regime, frequency divided by 2). The reference dashed horizontal lines are the reference frequencies: fres1 and fres2/2.

In the text

	Figure 4 Minimum value of the mouth pressure γ (green line) corresponding to a stable periodic solution (fundamental regime) with respect to the inharmonicity Inh (case Z1 > Z2; in practice Z1/Z2 = 3/2). Linear threshold γthr1 (blue dashed line).
In the text

	Figure 5 Signal (dimensionless acoustical pressure) obtained by simulation in time domain with an inharmonicity Inh = 0.04 and Z1 = 1.5Z2 (like in Fig. 3c). The dimensionless mouth pressure is printed in black, and is increasing linearly from γ = 0.43 (constant before t = 2 s) to γ = 0.50 (constant after t = 9 s).
In the text

	Figure 6 Signal (dimensionless acoustical pressure) obtained by simulation in time domain with an inharmonicity Inh = 0.04 and Z1 = 1.5Z2 (like in Fig. 3c). The dimensionless mouth pressure is printed in black, and is decreasing linearly from γ = 0.55 (constant before t = 2 s) to γ = 0.50 (constant after t = 9 s).
In the text

Figure 7 Three bifurcation diagrams of the two quasi-harmonic resonance air column (with Z1 < Z2, in practice Z2/Z1 = 3/2) with respect to the control parameter γ. From top to bottom: Inh = 0 (a), Inh = 0.015 (b) and Inh = 0.030 (c). Each case is described with two plots. Upper plot (a): Maximum of the absolute value of the periodic solution with respect to γ. The branch in red (respectively in blue, and in green) corresponds to the octave regime (respectively the fundamental regime associated to the standard Helmholtz motion, and the inverted Helmholtz motion fundamental regime). Lower plot (c): Frequency with respect to γ. The frequency branch in red (respectively in blue) corresponds to the fundamental frequency of the octave regime (respectively the fundamental regime). The frequency of the red branch has been divided by 2 for sake of clarity. The reference dashed horizontal lines are the reference frequencies: fres1 and fres2/2.

In the text

	Figure 8 Minimum value of the mouth pressure γ (blue line) corresponding to a stable periodic solution (fundamental regime) with respect to the inharmonicity Inh (case Z1 < Z2; in practice Z2/Z1 = 3/2). Linear threshold γthr1 (green dashed line).
In the text

	Figure 9 Signal (dimensionless acoustical pressure) obtained by simulation in time domain with an inharmonicity Inh = 0.015 and Z1 = Z2/1.5 (like in Fig. 7b). The dimensionless mouth pressure is printed in black, and is increasing linearly from γ = 0.38 (constant before t = 2 s) to γ = 0.45 (constant after t = 9 s).
In the text