Issue 
Acta Acust.
Volume 5, 2021



Article Number  33  
Number of page(s)  16  
Section  Musical Acoustics  
DOI  https://doi.org/10.1051/aacus/2021026  
Published online  09 August 2021 
Scientific Article
Multistability of saxophone oscillation regimes and its influence on sound production
^{1}
Aix Marseille Univ, CNRS, Centrale Marseille, LMA UMR7031, 13453 Marseille, France
^{2}
Laboratoire de Mécanique des Structures et des Systèmes Couplés, Conservatoire National des Arts et Métiers, 75003 Paris, France
^{*} Corresponding author: colinot.tom@gmail.com
Received:
26
August
2020
Accepted:
30
June
2021
The lowest fingerings of the saxophone can lead to several different regimes, depending on the musician’s control and the characteristics of the instrument. This is explored in this paper through a physical model of saxophone. The harmonic balance method shows that for many combinations of musician control parameters, several regimes are stable. Timedomain synthesis is used to show how different regimes can be selected through initial conditions and the initial evolution (rising time) of the blowing pressure, which is explained by studying the attraction basin of each stable regime. These considerations are then applied to study how the produced regimes are affected by properties of the resonator. The inharmonicity between the first two resonances is varied in order to find the value leading to the best suppression of unwanted overblowing. Overlooking multistability in this description can lead to biased conclusions. Results for all the lowest fingerings show that a slightly positive inharmonicity, close to that measured on a saxophone, leads to first register oscillations for the greatest range of control parameters. A perfect harmonicity (integer ratio between the first two resonances) decreases first register production, which adds nuance to one of Benade’s guidelines for understanding sound production. Thus, this study provides some a posteriori insight into empirical design choices relative to the saxophone.
Key words: Saxophone / Sound synthesis / Harmonic balance method / Numerical continuation / Attraction basins
© T. Colinot et al., Published by EDP Sciences, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
A classic endeavor in musical acoustics consists in the systematic study of sound production features of a musical instrument. Early studies use an artificial mouth to replace the musician (on the clarinet [1–4] or the bassoon [5]) in order to better describe and understand the physical phenomena at play during sound production. Later on, artificial mouths have been robotized to provide a complete mapping of the instrument’s behavior, aiming at understanding how the instrument must be acted on to produce different sounds [6, 7] or describing the influence of an acoustical parameter of the resonator on sound production [8]. The objective of this last study is shared by other works using a rather different approach to systematic description of the instrument’s behavior: using a physical model. Based on oscillation thresholds for instance [9], some conclusions can be drawn as to the acoustical characteristics facilitating the production of sound. Numerical resolution of the model’s equations also constitute a repeatable way to map the produced sound to the characteristics of the instruments, which has direct applications in instrument making [10, 11]. However, from a mathematical perspective, as nonlinear dynamical systems, wind instruments models often admit multiple solutions for a given set of parameters. The question of the stability of each of these solutions holds great importance when aiming to describe or predict the playability of an instrument based on its physical model. But some important questions remain unanswered, even for ideal cases where the stability or instability of each regime would be known. For instance, which regime is produced if two regimes are stable for the same control parameters combination?
In the case of such coexistence of stable solutions, denominated multistability hereafter, the convergence towards one or the other solution depends on the initial conditions. Indeed, each solution is associated with a region of attraction or attraction basin, defined as the region of the phase space where all initial conditions converge towards this solution [12, 13]. For instance, attraction basins are studied in walking models [14, 15], where the “walking” (periodic) regime almost always coexists with a stable equilibrium, corresponding to falling. In this case, describing attraction basins informs control strategies in robotics [16, 17]. Attraction basins are also studied for classic dynamical oscillators, such as Chua’s circuit [18], with experimental explorations of the attraction basins [19] as well as numerical investigations [20]. As strongly nonlinear selfoscillating systems capable of multiple oscillating regimes, wind instrument models are among the systems for which studying attraction basins can shed light on their rich behavior and help understand control strategies used by musicians. However, to our knowledge, no study on the attraction basins of musical instruments has been produced, although several studies explore their multistability. Experimental work on the clarinet [21] and a numerical study of several idealized woodwind resonators [22] illustrate in particular the hysteresis between regimes, which is a consequence of multistability. On the flute, continuation and synthesis have been used to investigate the hysteresis between regimes, notably depending on inharmonicity [23].
Describing the attraction basins and comparing their sizes is expected to give information on which regime is most likely produced, assuming some probabilistic repartition of the initial conditions in the phase space [24]. However, an exhaustive description is almost impossible for a complete model of instrument, where the phase space is of very large dimension. In such cases, attraction basins may be partially explored, based on a reduction of the phase space to one or two dimensions. For instance, the infinitedimensional phase space of a delayed system can be partially described along two dimensions [25]. In the case of musical instruments, a reduction of the phase space is proposed in this paper, based on knowledge of typical musical scenarios. Throughout this work, the case of a model of saxophone is considered, and two scenarios are studied: transition from another established limitcycle (scenario number 1), and first attack transient of a note, where the blowing pressure parameter goes from 0 to a certain final value (scenario number 2).
Section 2 presents the physical saxophone model and the two numerical methods used to solve its equations: the harmonic balance method and timedomain synthesis. Next, multistability is introduced by computing the bifurcation diagram with the harmonic balance method and continuation (asymptotic numerical method) and exhibiting hysteresis cycles using timedomain synthesis in Section 3 (control scenario number 1). Then, in Section 4, a simple testcase of scenario number two is presented to study sound production, where the blowing pressure increases from 0 to its final value over different durations. We show how this duration can influence the final regime in multistability regions, and explain these results by presenting the attraction basin of each regime. Section 5 demonstrates how the awareness of multistability can lead to a better description of the behavior of the model. Depending on the inharmonicity of the resonator, the size of the control parameter regions where each regime appears in synthesis is described, taking into account multistability. This provides an interpretation to the inharmonicity value measured on the saxophone by showing that it corresponds to an optimum in periodic regime production.
2 Numerical simulation framework
2.1 Saxophone model
The saxophone model used in this study is comprised of three main elements: a one degreeoffreedom oscillator representing the reed, a regularized nonlinear characteristic giving the flow through the reed channel, and a modal description of the measured impedance of the resonator. Similar models solved by timedomain synthesis (Sect. 2.2) are used in conjunction with analytical techniques to study the playing frequency [26] and spectrum [27] of clarinets, as well as their radiated power with a comparison to measurements [28]. The harmonic balance method (Sect. 2.3) can also be applied to this model to study its dynamic behavior, for instance to quantify the effect of neglecting reed contact [29].
Dimensionless [5, 30] acoustical Kirchhoff variables (p, u) are used in this work:
$$p=\frac{\widehat{p}}{{p}_{M}},\hspace{1em}u={Z}_{c}\frac{\widehat{u}}{{p}_{M}},$$(1)
where the hat notation indicates the variable in physical unit, p_{M} is the static pressure necessary to close the reed channel completely and Z_{c} is the characteristic input impedance of the resonator for plane waves. Similarly, the reed displacement from equilibrium is given in dimensionless form:
$$x=\frac{\widehat{x}}{H},$$(2)
where H is the distance between the reed and the mouthpiece lay at rest. With this formalism, the reed channel is closed when x ≤ −1. In this work, the only timevarying control parameter [4] is the dimensionless blowing pressure γ:
$$\gamma \left(t\right)=\frac{{p}_{m}\left(t\right)}{{p}_{M}},$$(3)
where p_{m} is the physical value of the pressure in the mouth of the musician. We leave all other control parameters constant in order to limit the dimensionality of the study. The values and names of the parameters are summarized in Table 1 and detailed below through the model description. Their values are drawn from [31] for the reed parameters q_{r} and ω_{r}, from [32] for the order of magnitude of the contact stiffness K_{c}.
Parameters of the numerical model: musician control parameters γ and ζ, reed parameters q_{r} and ω_{r}, contact parameter K_{c} and parameters inherent to the numerical implementation η, N_{m} and H.
2.1.1 The reed model
Following [32], the reed is modeled by a single degree of freedom oscillator including a nonlinear contact force accounting for the mouthpiece lay:
$$\frac{\stackrel{\u0308}{x}}{{\omega}_{r}^{2}}+\frac{{q}_{r}}{{\omega}_{r}}\stackrel{\u0307}{x}+x=p\gamma +{F}_{c}(x+1),$$(4)
where the two parameters of the reed are its angular eigenfrequency ω_{r} and its damping coefficient q_{r}, and the contact force is a function of the dimensionless reed opening x + 1 and is taken from [33],
$${F}_{c}(x+1)={K}_{c}\mathrm{min}(x+\mathrm{1,0}{)}^{2},$$(5)
where K_{c} = 100. Since x + 1 is the distance between the reed and the mouthpiece lay, F_{c} can be interpreted as a quadratic stiffness activated whenever the reed touches the lay. The ramp function min(x + 1, 0) is regularized using a parameter η = 10^{−3} to avoid nondifferentiability at x = −1 (reed closure),
$$\mathrm{min}(x+\mathrm{1,0})\simeq \frac{x+1\sqrt{(x+1{)}^{2}+\eta}}{2}.$$(6)
The regularization controlled by parameter η is necessary for the system to fit the quadratic formalism required by the implementation of the harmonic balance method and asymptotic numerical continuation in the MANLAB software, which produces the bifurcation diagrams of the present article. Although the parameter η is not necessary for the timedomain synthesis method to function, it is kept for comparison purposes.
2.1.2 The reed channel
The flow at the input of the resonator is deduced from Bernoulli’s law [34, 35] applied to the reed channel and turbulent mixing into the mouthpiece,
$$u=\zeta \mathrm{max}\left(x+\mathrm{1,0}\right)\mathrm{sign}\left(\gamma p\right)\sqrt{\left\gamma p\right},$$(7)
where ζ is the dimensionless control parameter accounting for reed opening at rest,
$$\zeta ={Z}_{c}\mathrm{wH}\sqrt{\frac{2}{\rho}},$$(8)
w being the effective width of the reed channel and ρ the density of the medium. We choose to ignore the flow due to the speed of the reed [26, 36] in the present model, as it only has a small effect on the playing frequency, which is not discussed here. The absolute value and ramp function in Equation (7) are regularized with the same parameter η as in Equation (6),
$$\gamma p\simeq \sqrt{(\gamma p{)}^{2}+\eta},$$(9)
$$\mathrm{max}(x+\mathrm{1,0})\simeq \frac{x+1+\sqrt{(x+1{)}^{2}+\eta}}{2}.$$(10)
2.1.3 The resonator
The input impedance is used to represent the resonator’s acoustical response. The dimensionless input impedance Z(ω) of a Buffet–Crampon Senzo alto saxophone is measured with the CTTM impedance sensor [37]. The saxophone is measured without its mouthpiece, placing the reference plane of the impedance measurement at the input crosssection of the crook. A cylindrical tube is added by transfer matrix method [38] in postprocessing to represent the mouthpiece, before using it in synthesis. The length of the cylinder is 60 mm and the radius is the same as the input radius of the crook, 6 mm. The total volume of the added cylinder approximately fits that of the missing cone apex, as per a classical academic approximation [39]. In order to use this input impedance for the two numerical synthesis methods presented above, it is decomposed into modes [26] so that,
$$Z\left(\omega \right)=\sum _{n=1}^{{N}_{m}}\mathrm{}\frac{{C}_{n}}{\mathrm{j\omega}{s}_{n}}+\frac{{C}_{n}^{\mathrm{*}}}{\mathrm{j\omega}{s}_{n}^{\mathrm{*}}},$$(11)
where C_{n} and s_{n} are the estimated complex modal residues and poles [40] and N_{m} is the number of modes retained in the simulation. In this paper N_{m} = 8 modes are used. This translates into the time domain by describing the pressure as a sum of complex modal components p_{n}, whose evolution depends on the modal coefficients, such that,
$${\stackrel{\u0307}{p}}_{n}\left(t\right){s}_{n}{p}_{n}\left(t\right)={C}_{n}u\left(t\right),\hspace{1em}\forall n\in [1,{N}_{m}],$$(12)
$$p\left(t\right)=2\sum _{n=1}^{{N}_{m}}\mathrm{}\mathrm{Re}\left({p}_{n}\right(t\left)\right).$$(13)
The flow u in (12) is given by (7). Figure 1 displays the measured impedance and the associated modal reconstruction according to Equation (11) for the D♯ fingering used throughout the rest of this article. The corresponding modal coefficients C_{n} and poles s_{n} are summarized in Table 2. Note that the choice of a modal formalism over a direct resolution of partial differential equations in the resonator is made here because of its lower computational cost and number of variables, which facilitates largescale numerical studies such as those presented in Section 5. Additionally, the modal formalism involves a limited number of parameters, which are directly tied to the acoustics of the resonator, instead of the complete geometrical description that a wave propagation model would require.
Figure 1 Measured input impedance (solid) and modal reconstruction (dashed) for the fingering D♯ of the alto saxophone. 
Modal parameters for the fingering D♯ of the alto saxophone.
2.2 Timedomain synthesis
Equations (4), (7) and (12) are discretized using finitedifference approximations for the timedomain derivatives, following a discretization scheme first applied to simple waveguides [41] and then adapted to a modal formalism [26]. The reader can find a detailed description of this discretization scheme in a recent document [42]. The sampling rate used in the simulation is F_{s} = 176 400 Hz, four times higher than the standard audio sampling rate. Such a high sampling rate is required, given the chosen finite difference scheme, to give precise results that match those obtained with the harmonic balance method.
As an illustrative result, Figure 2 shows an example of the synthesized pressure signal and its spectrogram. Note that the signal shown is a portion of the signal used in Figure 4, with control scenario number one. It corresponds to the first occurrence of the oscillations at a blowing pressure value γ ≃ 0.45. At this point, the system jumps from equilibrium to the first register and passes through fleeting second register and quasi periodic regimes. The spectrogram (Fig. 2b) shows the second register to be the octave (double the fundamental frequency) of the first register. The quasiperiodic portion of the signal displays amplitude variations, seen in the envelope of the signal (Fig. 2a) and on the odd harmonic components of the spectrogram. Quasiperiodic regimes are wellknown on saxophonelike instrument models, documented for instance in [2, 8, 43].
Figure 2 Timedomain synthesized pressure signal. (a) temporal envelope (black) and blowing pressure parameter γ (red). (b) Normalized spectrogram (dB) with regime names indicated (unstable ones between parentheses): equilibrium (Eq.), second register, quasiperiodic and first register. The signal is extracted from the same blowing pressure ramp as in Figure 4, between γ = 0.45 and γ = 0.51 (at first occurrence of oscillation). 
2.3 Harmonic balance and numerical continuation
The harmonic balance method is an analysis method particularly adapted to the study of musical instrument models [44], since it focuses on periodic solutions, which correspond to the produced notes. Assuming periodicity of the solution allows expanding all variables in Fourier series [45, 46] up to order H, such that the ith variable X_{i} expands to
$${X}_{i}\left(t\right)\simeq \sum _{h=H}^{H}\mathrm{}{X}_{i,h}\mathrm{exp}\left(\mathrm{jh}{\omega}_{0}t\right),$$(14)
where there are 2H + 1 complex Fourier coefficients X_{i,h} per variable, and ω_{0} is the fundamental angular frequency of the signal. In this work, the number of harmonics retained is H = 20. Applying the method to a differential system transforms it into an algebraic system of which the unknowns are the Fourier coefficients of the variables and the solution’s fundamental frequency, of the form,
$$R\left(\right\{{X}_{i,h}\},{\omega}_{0})=0.$$(15)
A numerical continuation method such as the asymptotic numerical method can then be applied to the resulting algebraic system [47, 48] to find how the solution changes for other constant values of a chosen control parameter, for instance the blowing pressure parameter γ. More precisely, knowing a solution to the system for a given value of γ, the continuation method finds a solution for a slightly higher or lower values of γ, and therefore progressively maps out the solutions of the system for a given range of blowing pressures. In this work, simulations were carried out using the MANLAB software (http://manlab.lma.cnrsmrs.fr/). This yields the value of the Fourier coefficients of the oscillating solutions along several values of a control parameter. The Fourier coefficients can then be used to reconstruct the timedomain solutions. This evolution can be summarized by a bifurcation diagram, which represents the variation of some descriptor, for instance the amplitude, of the solutions of the system with respect to the chosen control parameter. In addition, the stability of the solutions is determined using Floquet theory (for more details refer to [49–51]).
3 Multistability
This section presents the blowing pressure ranges where the model can produce each regime by studying their stability with the harmonic balance method. This result is summarized in the bifurcation diagram on which multistability zones appear as intervals where several regimes are stable. Signals are also synthesized with timedomain synthesis to exhibit how multistability leads to hysteresis. The correspondance between the two methods (the harmonic balance method and timedomain synthesis) is also checked.
3.1 Overlapping stability zones on the bifurcation diagram
The bifurcation diagram is computed for the (written) low D♯ fingering of an alto saxophone. The written D♯ produces the heard note F♯3 at 185 Hz as its first register regime. This intermediate fingering of the first register is chosen as test case because it exhibits both first and second register regimes, but no stable third register regime and few double twostep phenomena [52]. Figure 3 shows the L^{2}norm of the pressure signal:
$$\left\rightp{}_{2}=\frac{1}{T}\sqrt{{\int}_{0}^{T}\mathrm{}p(t{)}^{2}\mathrm{d}t},$$(16)
where T is the period of the signal, and identifies which regime each branch corresponds to. The blowing pressure parameter γ spans the interval between 0 and 2. This bifurcation diagram contains branches corresponding to the socalled equilibrium, where no sound is produced, for the lowest and highest γ values. The equilibrium at low γ corresponds to the musician not blowing hard enough into the instrument to obtain a sound, while at high γ equilibrium means that the reed channel remains closed due to the large pressure difference between mouth and mouthpiece. For intermediate γ values, the first and second register both appear. The first register is the fundamental pitch obtained with a given fingering, and the second register, sometimes referred to as overblowing, is pitched one octave higher than the first register. Even though the saxophone has an octave key facilitating the production of the second register, musicians know how to produce second register regimes without activating it. Note that contrary to the clarinet, the register key of the saxophone controls the two register holes, opening one or the other depending on all the pressed keys on the instrument. It is therefore not surprising that both regimes appear on the same fingering. In the present case, the γ interval between 0 and 2 contains all the studied limit cycles of the model, and at its bounds, only the equilibrium solution exists and is stable. The diagram in Figure 3 displays several zones of coexistence between stable regimes (i.e. multistability). These regions of coexistence are often bounded by bifurcation points, that mark qualitative changes in the oscillating regimes. In the present work, the system encounters several types of bifurcations, which we define succinctly in terms of the regime changes they correspond to. More formal definitions and characterizations of these bifurcations, notably in terms of critical values of the eigenvalues of the jacobian matrix of the system, can be found in [53, 54]. The Hopf bifurcation marks the emergence of an oscillating solution from equilibrium. The fold bifurcation corresponds to a stable and unstable solution branch colliding and disappearing, which can be better seen on the bifurcation diagrams as limit points in the solution branches. Neimark–Sacker bifurcations correspond to a periodic regime becoming unstable and being replaced by a quasiperiodic regime. A degenerate case of the Neimark–Sacker bifurcations is the perioddoubling bifurcation, where a periodic solution of halved frequency emerges from an oscillating solution. On saxophone models, period doubling bifurcations transform a second register regime into a first register. As is discussed in the next paragraph for the saxophone, Hopf and fold bifurcations often delimit coexistence between the equilibrium and an oscillating regime, while Neimark–Sacker and perioddoubling bifurcations mark the limits of multistability regions between two oscillating regimes.
Figure 3 Bifurcation diagram obtained with harmonic balance method and numerical continuation: L^{2}norm of the acoustical pressure depending on the blowing pressure parameter γ for the low written D♯ fingering of an alto saxophone. Thick lines: stable solution, thin lines: unstable solutions. Black: equilibrium, green: 1st register regimes, red: 2nd register regimes. Multistability zones are shaded: light yellow where 1st and 2nd register coexist, darker gray for equilibrium and 1st register. Blue circles specify the location of bifurcations. Vertical black lines correspond to those in Figure 6c, and (from left to right) to phase diagrams 7a–7c. 
Starting from low blowing pressure values, the first coexistence zone appears between the first register and the equilibrium. It is delimited by the fold bifurcation F1 of the first register around γ = 0.3 and the inverse Hopf bifurcation H1 at γ = 0.4 where the equilibrium becomes unstable. The second coexistence zone is between first and second register, in the interval where the second register is stable between the Neimark–Sacker bifurcation NS1 and the perioddoubling bifurcation PD1 (respectively at γ = 0.66 and γ = 0.79). The Neimark–Sacker bifurcation NS1 mark the destabilization of the second register and the emergence of a quasiperiodic regime (not represented here), sometimes called multiphonics by musicians. The next coexistence zone occurs in the interval between the two perioddoubling bifurcations PD1 and PD2 on the second register branch, where a stable double twostep solution [52] emerges. This coexistence zone is not shaded on the figure, as it could represent less of a musical issue, since double twostep regimes have roughly the same frequency as standard first register regimes. The fourth coexistence zone is more complicated: it starts between first and second register at the perioddoubling bifurcation PD2, and then the equilibrium also becomes stable at the Hopf bifurcation H4. The limit of the last coexistence zone is made by the two fold bifurcations F2 and F3 where the first and second register solutions cease to exist.
The diagram in Figure 3 shows that coexistence zones between stable regimes span most of the range in γ where oscillating solutions exist, including arguably crucial γ values like the lowest for which an oscillating regime exists. Multistability is not an isolated phenomenon, but rather corresponds to the general situation, at least for this fingering.
3.2 Timedomain synthesis with blowing pressure ramps (control scenario number one)
Once multistability zones are identified, timedomain synthesis can be used to exhibit their role when playing the instrument. One of the main phenomena multistability entails is hysteresis: for several values of the blowing pressure, a different regime is produced depending on whether the blowing pressure is increasing or decreasing. Various multistable regimes are exhibited using this method in [22] on woodwind models. Figure 4 shows the hysteresis cycles obtained by using ramps of γ: in control scenario number one, the parameter γ is progressively increased from 0 to 2 and then decreased back to 0. Each one of the increasing and decreasing phases of the synthesis has a duration of 60 s. This duration was chosen after several trials, sufficiently long to let stable regimes establish while keeping a γ slope steep enough to limit dynamical bifurcation delays [55].
Figure 4 L^{2}norm of the timedomain synthesis signal (dark blue line) for control scenario number one, a ramp of the blowing pressure parameter γ from 0 to 2 and down, overlaid with the bifurcation diagram of Figure 3 (thick stable solutions and thin unstable solutions, of the first register in green and second register in red). Solid dark blue line indicates the first part of the signal (increasing γ) and dashed dark blue line the second (decreasing γ). 
Figure 4 shows that the synthesized signal starts from γ = 0 at equilibrium, its L^{2}norm being zero. Then, at Hopf bifurcation H1, the equilibrium becomes unstable, which causes the system to start oscillating. At this point, the synthesis goes through a transient represented in Figure 2 (between t = 14.2 s and t = 14.8 s), shortly passing by unstable second register and quasiperiodic regimes before reaching the first register. Once the first register is established, the branch is followed all the way to extinction (around γ = 1.75 on Fig. 4), because the first register does not become unstable until the fold bifurcation F3. At this point, the system returns to equilibrium until the highest value of the ramp, γ = 2. The blowing pressure γ then starts descreasing, and the system stays at equilibrium until Hopf bifurcation H3 is reached, for γ ≃ 1.05. There, the system jumps to the stable second register regime. The second register branch is followed until the perioddoubling bifurcation PD2, where the system briefly follows the double twostep branch. The double twostep branch appears on the L^{2}norm very close to the second register branch, and looks like a small bump. The two branches cross again at perioddoubling bifurcation point PD1. Then, something rather surprising occurs: the L^{2}norm of the timesynthesized signal seems to follow the second register branch of the bifurcation diagram further than Neimarck–Sacker bifurcation NS1, although the second register branch has become unstable. This is because the quasiperiodic regime emerging at NS1 is actually a stable attractor, and the associated L^{2} norm happens to be close to that of the second register. Branches of quasiperiodic regimes have not been computed so as not to clutter Figure 3, but note that it is possible with harmonic balance method and MANLAB [48]. When the quasiperiodic regime becomes unstable (around γ = 0.55) the system jumps back onto the first register branch, which is followed until fold bifurcation F1, beyond which the stable equilibrium is the unique solution of the model. The path described precedently is highly hysteretic: the sequence of regimes produced for increasing and decreasing γ are very different. Actually, the two paths only coincide in three regions: the lowest and highest γ intervals, for which only the equilibrium is stable, and a very small region around γ = 0.5 where only the first register is stable.
The hysteresis phenomenon observed here in timedomain synthesis can be interpreted as the first step in attraction basin description: once a certain stable regime is reached, it is followed until extinction or loss of stability, even when other regimes are simultaneously stable. This confirms that a stable periodic regime is part of its own attraction basin. Reconstructing stable parts of the bifurcation diagram using timedomain synthesis and comparing them to those obtained using the harmonic balance method also provides validation for the numerical discretization scheme in this context. Here, it shows that timedomain synthesized signals are not perturbed by numerical artifact due to time discretization and can be trusted to describe properties of the model. Note that even subtle details of the behavior such as the tiny branch of double twostep solution between period doubling PD1 and PD2 are found by the timedomain synthesis.
This exploration of the blowing pressure space using a long ramp is very useful to exhibit the hysteresis phenomenon, as well as test the coherence between the two synthesis methods. However, this kind of sound is extremely artificial and far from anything a musician would use in everyday practice (provided it is even possible for a musician to produce it). Therefore, we frame the conditions of the rest of the study so that they can be interpreted in terms of selection of one regime over another.
4 Regime emergence in a multistable context
It is very likely that musicians learn to select between coexisting stable regimes, adjusting their control so that the established regime in a multistability region is the one they desire. This idea provides the layout for control scenario number two: we study the effect of a parametrized transient control or initial conditions on the established steadystate regime that follows when the control is constant. To provide another way to interpret the selection process between multistable regime, a third control scenario is used, where only initial conditions are varied and all control parameters are constant.
4.1 Control scenario number two: increasing blowing pressure
One way to study the attraction basins more thoroughly is to run many simulations with initial conditions spanning the whole phase space. However, since the considered model has a 2N_{m} + 2 dimensional phase space, a complete exploration is not possible. Moreover, many of the possible initial conditions are unlikely to be created by the musician. More interesting is the exploration of the regions of the phase space that are crossed by the system when a given control pattern is varied. Here, we focus on a monotonic increase of the blowing pressure γ at the attack: without using the tongue, the player starts blowing progressively harder into the instrument. Such a scenario was proposed in [56]. Note that instrumented mouthpiece measurements performed on the saxophone, such as those presented in [57], often show a different profile including a pressure overshoot before the apparition of the oscillations. However, in the present study, we omit that overshoot so that the control scenario is entirely defined by a single parameter. In control scenario number two, the blowing pressure starts from 0 and rises up until stabilizing at a certain value γ_{f}, during a certain time determined by the parameter τ_{g}. The temporal variation of γ is given by the expression,
$$\gamma \left(t\right)=\frac{{\gamma}_{f}}{2}\left(1+\mathrm{tanh}\left(\frac{t5{\tau}_{g}}{{\tau}_{g}}\right)\right),$$(17)
which is differentiable infinitely many times. Figure 5 displays four examples of such transients.
Figure 5 Four examples of blowing pressure variations in control scenario number two, as described by Equation (17), for γ_{f} = 0.4 and 0.6 and τ_{g} = ${\tau}_{g}^{\left(1\right)}$ = 10 ms and ${\tau}_{g}^{\left(2\right)}$ = 20 ms. 
Other envelopes (sigmoid, sine branch) were tested, only causing small quantitative changes to the results. Figure 6 shows which established regimes appear in timedomain synthesis depending on τ_{g}, for final values γ_{f} belonging to the multistability zones described in Figure 3. Each dot on the figure represents the type of established regime after five seconds of timedomain synthesis. This synthesis duration was chosen to be sufficiently long so that the transient is completed and the established regimes can be observed. Regime types are estimated using an energybased criterion for equilibrium (if the energy of the pressure signals in the last ten periods is less than that of the first ten, the regime is classified as equilibrium) and a fundamental frequency estimator to determine the first and second register. To detect quasiperiodic regimes, an attribute that can be observed in Figure 2 is used: the fact that the amplitudes of the harmonics of a quasiperiodic signal vary temporally. For the classification, if the mean variance of the harmonic amplitudes is more than a certain threshold (here set to 10^{−6}), then the regime is considered quasiperiodic.
Figure 6 Classification of the steadystate regimes produced depending on the blowing pressure transient parameters: final value γ_{f} and characteristic time τ_{g}. Multistability zones (deduced from Fig. 3): (a) equilibrium ○ and first register (b) first register and second register (c) all three regime types. Blue triangles indicate quasiperiodic regimes. The horizontal line on graph (b) shows τ_{g}= 1/f_{1}. The vertical lines highlight the γ values of phase diagrams in Figures 7 and 8. 
Figure 6a focuses on the first multistability region (highlighted in gray in Fig. 3), near the first Hopf bifurcation H1. The two stable regimes in this region are the equilibrium (○) and the first register (). For final values γ_{f} between 0.38 and 0.4, the system can converge to both regime depending on the characteristic rising time τ_{γ}. It is interesting to note that equilibrium is reached for the longest rising times, i.e. the slowest γ variation, whereas the oscillating regime is reached for the shortest rising times. This is understandable as a quick γ increase tends to drive the system away from equilibrium, and therefore possibly out of its attraction basin. Note that some of the oscillating regimes near the limit, for the longest attack times, are classified as quasiperiodic. This is due to the transient being extremely long in this particular region: the steadystate regime is not yet established at the end of the synthesized signals. This classification particularity, which could be seen as an error, is not corrected because, from a musician’s perspective, a regime still varying five seconds after the start of the attack will arguably not be considered periodic.
The second zone of multistability is explored in Figure 6b. The first and second register are separated by some stable quasiperiodic regimes (). This is the same quasiperiodic regime that appears in timedomain synthesis in Figure 2, which overlays the unstable portion of the second register branch. There is a particular range of characteristic time τ_{g} that seems to produce the first register for a larger range of γ_{f}. The corresponding values of τ_{g} are close to the period of the first register, represented by the horizontal line in 5 (b).
In the last multistability (0.9 ≤ γ ≤ 1.25), three regimes may be stable for the same parameter values, as shown on Figure 3. However, Figure 6 reveals that there is no γ_{f} region where all three are produced. This can be explained by analyzing the attraction basins (see Sect. 4.2, Fig. 8).
4.2 Control scenario number three: varying initial conditions
The results concerning the influence of the blowing pressure parameter can be better understood by examining the region of the phase space leading to each regime. A point in the phase space represents the current state of the system, meaning the value of the state variables and their derivatives. Since the system is deterministic, a given point in the phase space will always lead to the same stable established regime. Therefore regions of the phase space can be associated with each regime. These regions are called attraction basins. The attraction basins are only rigorously defined in a context where all control parameters are constant. However, to interpret the results obtained with a control parameter transient such as control scenario number two, it is useful to observe the layout of attraction basins obtained for the final blowing pressure value γ_{f}. Specifically, at the beginning of control scenario number two, the blowing pressure is subject to fast variations, making a direct approach based on attraction basins illdefined. However, after the transient, the blowing pressure stabilizes around its final value γ_{f}. To elucidate the behavior of the system at this moment, Section 4.2 performs a systematical analysis of its convergence with constant control parameters, in which case the attraction basins representations are relevant. Therefore, a third control scenario is devised, where the control parameters are kept strictly constant and only the initial values of certain state variables of the system are modified.
Because the phase space is of dimension 2N_{m} + 2 (all modal components and their derivatives, plus reed position and speed), it is necessary to choose a projection to represent the attraction basins. After some trials, a projection of the phase space on the two first modal components (see Eq. (12)) and the derivative of the second one, (p_{1}, p_{2}, ${\stackrel{\u0307}{p}}_{2}$), was chosen as a threedimensional projection. These variables were chosen not only because of their physical or mathematical meaning, as they relate respectively to the first and second register, but also because they allow for the clearest visual separation of the limit cycles and attraction basins that could be obtained by the authors. To estimate the attraction basins, timedomain synthesis is launched with initial conditions spanning the projected phase space and constant control parameters (control scenario number three) . A total of 256 initial conditions are scattered in a Latin hypercube sampling into a rectangular parallelepiped such that,
$${p}_{1}^{I}\in [\mathrm{0.2,0.2}],\hspace{1em}{p}_{2}^{I}\in [\mathrm{2,2}],\hspace{1em}{\stackrel{\u0307}{p}}_{2}^{I}\in [\mathrm{707,707}].$$(18)
These bounds should be understood with respect to the amplitude of the limit cycle along each dimension (that can be seen in Figs. 7 and 8). They were chosen so that whenever a regime is stable, it is obtained in synthesis at least once. All the other modal pressure components and their derivatives, as well as the reed speed $\stackrel{\u0307}{x}$, are initially zero. So that there is no discontinuity when starting the synthesis, the initial values of the variables p, then x and u are computed accordingly through Equations (13), (4) and (7). All control parameters, including the blowing pressure γ, are kept constant during these simulations: recall that the attraction basin analysis is only strictly valid for constant values of the control parameters. Therefore, the deductions that are made using the projected attraction basins (control scenario number three) on the results of control scenario number two are subject to caution.
Figure 7 Projection of the attraction basins obtained with control scenario number three. Large dots are initial conditions, small dots are points the synthesis goes through. The dots’ colors indicate the final regime they lead to (black: equilibrium, green: first register, red: second register). Blue lines represent the limit cycles. The small one on the inside is the first register and the one on the outside is the second register. The blowing pressure γ is (a) 1, (b) 1.1, (c) 1.2 (highlighted in Figs. 3 and 6). 
Figure 8 Attraction basins (dots) and limit cycles (lines) in a 3D projection of the phase space obtained with control scenario number three for different blowing pressures γ. Black: equilibrium, green: first register, red: second register, blue: quasiperiodic. 
Figure 7 shows these initial conditions in the plane (p_{2}, ${\stackrel{\u0307}{p}}_{2}$), associated with the regimes they lead to, for three values of the blowing pressure parameter γ. All the phase space points the system passes through during the transient are also part of the attraction basin, so they are represented in the figure as well. The last period of synthesized signals is plotted as the limit cycles. The three values of γ are chosen near the last Hopf bifurcation (H3 and H4 in Fig. 3), where the equilibrium becomes stable again. Graph 6 (a) is computed at γ = 1, before the Hopf bifurcations, so only register one and two are stable. Therefore, none of the initial conditions lead to equilibrium (no black dots on the figure). Each attraction basin is located around the corresponding limit cycle. The attraction basins seem to overlap, but it is merely an effect of the projection of the phase space. Graph 6 (b) corresponds to γ = 1.1, right above the Hopf bifurcations H3 and H4. One can see that some initial conditions located near the origin now lead to equilibrium. This plot can seem surprising when compared to Figure 6c, which only shows equilibrium and register one for this value of γ_{f} using the control scenario number two, although the attraction basin of the second register (red points in Fig. 7b) seems larger than that of the other regimes. This is due to two effects. First, the projection spreads the second register but shrinks the first register. Secondly, control scenario number two starts from the origin of the phase space: it is ensured by setting all variables and their derivatives at 0 at the start of the synthesis. Therefore, attraction basins surrounding the origin, such as that of the first register, are more likely to be entered. Here we see that considering only the size of the attraction basins and ignoring their position provides only little information. This data needs to be interpreted in terms of musician action, notably pertaining to the ability of the musician to actually reach the region of the phase space. Here, the fact that the attraction basin is vast is effectively countered by the fact that it is far from the origin, making it hard to reach using the studied control scenario number two. Graph 6 (c) represents the same results for γ = 1.2, where Figure 6c showed more occurences of the equilibrium than for γ = 1.1. This is explained by the attraction basin of the equilibrium being larger – there are many more black dots on Figure 7c than on (b). Notice that the attraction basin of the first register expands (also many more green dots than on Fig. 7b), while the attraction basin of the second register (red dots) shrinks. This process (not represented here) continues until the second register ceases to be stable at fold bifurcation F3 (see Fig. 3).
Figure 8 shows the attraction basins and limit cycles, in a threedimensional projection of the phase space (p_{1}, p_{2}, ${\stackrel{\u0307}{p}}_{2}$), at particular values of γ highlighted in Figure 6. Graphs 2 (a), (b), (d) and (e) should be read as further information on the regime map 5 (b), at the beginning of the multistability zone between first and second registers.
Graph 2 (a) corresponds to γ = 0.6, and confirms that the first register is the only stable regime: it was the only one to appear in the regime map 5 (b) (for γ_{f} = 0.6). Then, a quasiperiodic attractor appears in Graph 2 (b), for γ = 0.63. Although the associated attraction basin seems smaller than that of the first register, it seems to almost surround the origin of the phase diagram. It was observed in synthesis that the transient of control scenario number two does not send the system very far from the origin in the phase space – compared, for instance, with the size of the limit cycles of register one and two. This is linked to it necessarily starting from the origin of the phase space. The fact that control scenario number two tends to lead the system to phase space points near the origin explains why regime map 5 (b) displays more quasiperiodic regimes than first register. A similar interpretation can be formulated with regards to graphs 2 (d) and (e), for γ = 0.645 and γ = 0.72 respectively. For these values of γ on Figure 6b, there is more second register than first register. On graphs 2 (d) and (e), it can be seen that although the size of the first register’s attraction basin seems comparable to that of the second register, the latter clearly holds a central position around the origin of the phase space. When the second register attraction basin grows in Graph 2 (e), this translates to the disappearance of the first register from the regime map 5 (b). Note that Graph 2 (e) confirms that the first register is still stable, as announced by the harmonic balance method in Figure 3. Graph 2 (c) γ = 0.9 illustrates a slightly different explanation of a similar case of only second register appearing in the regime map 5 (c): in this case, the attraction basin of the first register is just too small, it only makes up for a few points in Figure 8c. Graph 2 (f) (γ = 1.25) is comparable to (d) in that many regimes are stable, but only the one with the most central attraction basin appears in the corresponding regime map on Figure 6c. In this case this regime is the equilibrium, whose attraction basin in Graph 2 (f) surrounds the origin, although it appears smaller than the others. To complete the study, the full evolution sequence of the attraction basins obtained with control scenario number three can be found as an animation in multimedia file Supp1.mp4. The authors suggest frequently pausing the animation to observe precisely how the attraction basins develop along multistability zones.
5 Effect of the resonator’s inharmonicity on regime production
The concept of multistability and the attraction basins are presented and explored here because they seem to be a very important part of the observed behavior of a saxophone model. The present section offers a succinct description of the behavior of the model, applied to an instrument design problematic. It also illustrates how ignoring multistability can affect the description of the behavior of a model.
Before using the analysis of a woodwind physical model to develop new instruments, it can be very informative to apply it to existing instruments, in the idea of a reverse engineering procedure. If the analysis method can explain a posteriori some design choices made on instruments with satisfying sound production characteristics, then it might help guide further innovative design choices in the right direction. In the present case, the produced regimes are studied for the seven lowest first register fingerings, and one acoustical parameter is varied artificially: the inharmonicity between first and second resonance. The original data corresponds to measured impedance for the corresponding fingerings of a Buffet–Crampon Senzo professional alto saxophone. According to the socalled Bouasse–Benade prescription [58–60], near perfect inharmonicity between the resonances is cited as a condition for good playability of the instrument. This prescription is also discussed in recent studies [9, 61]. On the saxophone, experimental studies using an artificial mouth have shown that varying inharmonicity greatly affects regime production [2, 8]. In this work we define inharmonicity as the ratio between the second and first resonance f_{2}/f_{1}, or $\mathfrak{I}\left({s}_{2}\right)/\mathfrak{I}\left({s}_{1}\right)$ in terms of the parameters of Equation (11). On a saxophone this ratio is close to two. Many definitions of the inharmonicity can be devised, possibly taking into account more resonances. The present definition of the harmonicity has the advantage of being very short, and easy to modify in the modal formalism by adjusting only one modal frequency.
Specifically for the purpose of the following study, optimal regime production conditions are defined crudely in terms of how often each regime appears in synthesis. This definition differs from that employed in [9]. On the lowest fingerings, in this study, we simply consider that optimizing regime production means maximizing the appearance of the first register while minimizing that of the second register and quasiperiodic regimes. Indeed, one of the challenges many beginner saxophone players face on the lowest fingerings is controlling the instrument so that the first register can be produced, and not another regime. Quasiperiodic regimes are largely considered undesirable in common musical practice, however they are a common issue on the lowest fingerings of the saxophone.
5.1 Regime production regions
Expanding on the idea in Figure 6, one can study the produced regime across the twodimensional parameter space (γ_{f}, ζ), while still varying the characteristic time τ_{g} of control scenario number two. Figure 9 shows the classification of obtained regimes for several combinations (γ_{f}, ζ) and several characteristic times. Note that for a musician using their lower lip to control the instrument, it is difficult to control the reed opening parameter ζ without also varying the properties of the reed ω_{r} and q_{r}. Leaving the reed parameters constant in this study amounts to a simplification of the musician control. For readability reasons, the resolution of the cartography presented here is rather coarse: only eight values of γ_{f} and ζ and three characteristic times τ_{g}, for a total of 192 synthesized signals. The range in parameter ζ is inspired by the range measured in [8] for artificial mouth experiments, using the method proposed in [36]. This method relies on measuring the flow rate as a function of the mouth pressure, and estimating ζ based on the maximum flow using the nonlinear characteristic of Equation (7). As a case study, two maps computed with different inharmonicity values are presented on Figures 9a and 9b, so that they can be compared. In the modal formalism, the inharmonicity is changed very simply by modifying the value of the second modal frequency. Note that changing the length of the saxophone mouthpiece constitutes another method to vary the inharmonicity. However, it modifies all the modal frequencies simultaneously, which makes some interpretations less robust. Therefore, only the results obtained by varying the second modal frequency are presented here. Similar results can be obtained by varying the length of the mouthpiece. Two typical values of inharmonicity are chosen: one that could be called null, f_{2}/f_{1} = 2; and the value measured on the saxophone which is slightly higher, f_{2}/f_{1} = 2.065.
Figure 9 Classification of the regimes produced (○ equilibrium, first register, second register, quasiperiodic) depending on control parameters: γ_{f} and ζ. Each rectangle corresponds to a couple (γ_{f}, ζ) and the points inside indicate the regime for each characteristic time τ_{g} (bottom 0.1 ms, middle 3 ms and top 100 ms). One rectangle on graph (b) is annotated as an example. Graphs correspond to two inharmonicities for low written D♯ fingering (a) f_{2} = 2.065f_{1}, the value measured on a real saxophone and (b) f_{2} = 2f_{1}. 
Focusing on Figure 9a, several features can be described, and recognized from the situations explored in Section 4 with a fixed ζ. Coexistence regions can be noticed on most of the map, with a given (γ_{f}, ζ) couple leading to different regimes depending on the characteristic time. This further demonstrates that multistability is a very common phenomenon across the control parameter space in woodwind models. A particular case of coexistence occurs on all regime maps, where long attack times lead to the system remaining at equilibrium, while fast attacks can trigger oscillations. This is the same phenomenon as in Figure 6a. These phenomena are located, as for Figure 6a, at the boundaries between equilibrium and oscillation regimes, which correspond to the Hopf bifurcations of the model (where the equilibrium becomes unstable). On the regime maps, points containing both equilibrium regimes and oscillating regimes are seen on the vertical threshold near γ = 0.4, the horizontal threshold near ζ = 0.2 and the extinction threshold around γ = 1.1.
Coexistence situations similar to Figure 6b can also be seen on Figure 9a, for example at ζ = 1.2 and γ_{f} ≃ 0.8, where the short and long characteristic times lead to the second register, while the medium time leads to a first register. This possibly indicates that the second register attraction basin almost surrounds the origin, as seen on Figure 8, in all the multistability zones between first and second register for that fingering.
Figure 9b (f_{2}/f_{1} = 2) displays a lot more second register () than Figure 9a (f_{2}/f_{1} = 2.065). Contrary to what could be expected, a null inharmonicity, where the second resonance frequency is twice the first, does not lead to more first register production. Since an exact integer ratio between resonances does not facilitate the production of the first register, one can ask if the model shows a particular value of inharmonicity which favors the production of first register.
5.2 Rate of produced regimes: influence of the rise time on global regime production
To study the question of the inharmonicity favoring first register production, regime maps are computed for all the fingerings of the saxophone that should produce first register regimes, meaning those where the register holes are closed. For each fingering, the second modal frequency f_{2} is varied from 1.96f_{1} to 2.15f_{1}, by steps of 0.01f_{1}. A regime map containing N_{p} = 192 points, as for Figure 9, is then computed for each value of inharmonicity. The produced regimes are counted for the whole map, and a rate is computed for each of them with respect to the total number of oscillating regimes as,
$${R}_{i}=\frac{{N}_{p,i}}{{N}_{p,o}},$$(19)
where regime i can either be first register, second register or quasiperiodic regimes, N_{p,i} is the number of points corresponding to regime i in the regime map and N_{p,o} is the total number of points corresponding to any oscillating regime (i.e. all regimes but equilibrium). Note that this description ignores nonoscillating regimes.
Figure 10 depicts the rate of each oscillating regime produced depending on inharmonicity for the lowest fingering of the first register, the written Bb which produces the heard note C3 at 131 Hz. On this figure the regimes are counted separately for each characteristic time τ_{g} before being summed for the whole map to produce an averaged rate. On the averaged rate, optimal points are highlighted by triangles, corresponding to the maximum of first register and minimum of second register produced, respectively. Both appear for values slightly above f_{2}/f_{1} = 2. The inharmonicity values maximizing first register production do not correspond to exactly harmonic resonances, but a second resonance slightly higher than the octave of the first. The proportion of quasiperiodic regimes is also displayed on the figure. Note that inharmonicity values around two lead to less quasiperiodic regimes, which is corroborated by existing results [8, 43]. On Figure 10, it can be seen that the region of low rate of production of the second register for f_{2}/f_{1} between 2.05 and 2.09 almost coincides with a region of high rate of production of quasiperiodic regimes, for f_{2}/f_{1} between 2.06 and 2.1. In this case, minimizing second register production seems to favor the production of quasiperiodic regimes.
Figure 10 Rate of produced regimes (Eq. (19)) for written low B♯ fingering. Green: first register, red: second register, blue: quasiperiodic. Linestyles indicate the characteristic time. Dotted: τ_{g} = 0.1 ms, dashdot: τ_{g} = 3.2 ms, dashed: τ_{g} = 100 ms, solid: averaged rate. An upward triangle marks the maximum first register averaged rate, a downward triangle marks the minimum second register rate. 
Figure 10 shows that ignoring multistability, by using only one value of characteristic time, could have lead to conclusions similar to those drawn with production rates averaged over several characteristic times. Indeed, all the characteristic times yield qualitatively similar production rates. However, this is not always the case, and Figure 11 shows two examples where neglecting multistability and studying only one characteristic time can lead to biased conclusions.
Figure 11 Rate of first register regime (Eq. (19)) produced for (a) written low D♯ fingering and (b) written low D fingering. Linestyles indicate the characteristic time. Dotted: τ_{g} = 0.1 ms, dashdot: τ_{g} = 3.2 ms, dashed: τ_{g} = 100 ms, solid: averaged rate. An upward triangle marks the maximums of first register rates. 
For the D♯ fingering (Fig. 11a), one can see that depending on the chosen increase duration τ_{g} the production ratio varies greatly (from 20% to 60%). If any quantitative interpretation is to be expected from these results, it can be changed dramatically depending on the chosen attack time. Notice that the lowest rate of first register corresponds to the longest attack time τ_{g} = 100 ms. Figure 9 shows that the first register regimes are produced on the edges of the zone of oscillation, in a multistability region between the first register and the equilibrium. The longest attack times in this region tend not to lead to a first register, but instead to an equilibrium due to its attraction basin surrounding the origin (see for instance Fig. 7c or multimedia file Supp1.mp4). Figure 11b shows the results for the written D fingering. This case exhibits an outlier: the shortest attack time yields a optimal inharmonicity value of 2.08, whereas the others point to 2.04. In this case, considering several attack times is a way to smooth out outliers due to a particular value of the attack time.
Note that Graph 9 (b) can also be subject to an interesting interpretation in terms of musician control strategies: the fact that certain attack time values seem to markedly decrease the rate of production of a certain regime could be used by the musician to avoid producing it.
5.3 Inharmonicity of the saxophone
In this section, the optimal inharmonicity in terms of regime production is studied for the seven lowest fingerings of the instrument. Higher fingerings are not represented because they add no relevant information: first register regime production rates are close to 100% for all the studied inharmonicity values. This corresponds to the saxophonists’ experience that the high notes of the instrument’s first register are often easier to produce than the low notes, and to the fact that the first impedance peak is much higher than the others on the high fingerings [62]. The optimums are compared with the inharmonicity value measured on the saxophone on which the model is based. Figure 12 summarizes the production ratios for all the fingerings. The optimal inharmonicity seems to vary across the fingerings. It is always greater than two: null inharmonicity does not favor first register production on the low fingerings of the saxophone. The two optimums are close to the measured inharmonicity. Additionally, the trend is respected, with optimal and measured harmonicities increasing for higher fingerings. Note that the optimum for the E fingering is very far from the measured inharmonicity, but the production ratios are almost constant. Overall, the simulation shows that the most first register and least second register is produced by the model for values of inharmonicity near those measured on the saxophone. This result sheds some light on the empirical choice of the acoustical properties on the saxophone. Indeed, if the inharmonicity was far from the values observed on saxophones, the model predicts more second register would be produced. This effect is arguably undesirable. However, this choice comes with a compromise, as it also favors quasiperiodic regime production (see Fig. 10), which are a known issue in low saxophone fingerings.
Figure 12 Rate of produced regimes for the lowest fingerings of the alto saxophone (written pitch). Green: first register, red: second register. An upward triangle marks the maximum first register averaged rate, a downward triangle marks the minimum second register rate. Vertical lines mark the measured inharmonicity values. 
6 Conclusion
In the studied saxophone model, stable regimes coexist throughout large regions in the musician control parameter space. Thus, even an exhaustive description of the stability or instability of each oscillating regime in the control parameter space is only an incomplete answer, as it doesn’t suffice to predict which regime emerges in these multistability zones. In particular, in the event that the attraction basin of a stable regime proves too small or unattainable by the musician, this regime should be considered unreachable, almost like an unstable regime. This work demonstrates that saxophone regime stability is a nuanced topic.
An adapted combination of two radically different numerical methods provides a more complete description of the model’s behavior: the stability study performed by harmonic balance method is completed by timedomain synthesis, which quantifies multistability by outlining the attraction basin of each regime. A varying control scenario explores which regime is produced in multistability zones, with the advantage that it can be tied to plausible musician actions. However, a varying control scenario only provides a very partial view of the attraction basins, and its results deserve to be explicited by representing the attraction basins in the phase space using different initial conditions. Dedicated experimental work, out of the scope of this paper, could help design more realistic control scenarios.
Accounting for multistability, the study of synthesized regimes may explain an acoustical choice made by instruments makers: the inharmonicity of the saxophone. An integer ratio between the first and second resonance frequencies does not favor the production of first register. Note that this result adds nuance one of Benade’s guidelines stating that an oscillation is favored if the impedance is large at its fundamental and its harmonic frequencies [59]. This work shows that competition between registers also comes into play, depending on more than solely the impedance magnitude at the playing frequency and its harmonics. Instead, an integer ratio between the first and second resonance frequencies tends to favor the production of second register, which is arguably undesirable for a first register fingering. Carefully tuned inharmonic resonances, where the second frequency is higher than twice the first, can lead to more first register production. The optimal inharmonicity value found on the model is close to harmonicities measured on saxophone resonators. This result provides an a posteriori interpretation of the acoustical characteristics of the saxophone, as chosen empirically by instrument makers, as the acoustical characteristic leading to easier production of the first register. Such results are among the first steps towards applying numerical simulations as predictive tools to estimate playability in instrument design.
Supplementary materials
Supplementary material is available at https://actaacustica.edpsciences.org/10.1051/aacus/2021026/olm
Supp1.mp4. Access here
Conflict of interest
Author declared no conflict of interests.
Acknowledgments
The authors would like to thank Erik Petersen for correcting the English. This work has been carried out in the framework of the Labex MEC (ANR10LABEX0092) and of the A*MIDEX project (ANR11IDEX000102), funded by the Investissements d’Avenir French Government program managed by the French National Research Agency (ANR). This study has been supported by the French ANR LabCom LIAMFI (ANR16LCV200701).
References
 J. Backus: Vibrations of the reed and the air column in the clarinet. The Journal of the Acoustical Society of America 33 (1961) 806–809. [CrossRef] [Google Scholar]
 J.P. Dalmont, B. Gazengel, J. Gilbert, J. Kergomard: Some aspects of tuning and clean intonation in reed instruments. Applied Acoustics 46 (1995) 19–60. [CrossRef] [Google Scholar]
 C. McGinnis, C. Gallagher: The mode of vibration of a clarinet reed. The Journal of the Acoustical Society of America 12 (1941) 529–531. [CrossRef] [Google Scholar]
 T.A. Wilson, G.S. Beavers: Operating modes of the clarinet. The Journal of the Acoustical Society of America 56 (1974) 653–658. [CrossRef] [Google Scholar]
 A.Y. Gokhshtein: Selfvibration of finite amplitude in a tube with a reed, in Soviet Physics Doklady, MAIK Nauka/Interperiodica and Springer Science+Business Media (Russia). 1979, Vol. 24, p. 739. [Google Scholar]
 T. Helie, N. Lopes, R. Causse: Openloop control of a robotized artificial mouth for brass instruments. The Journal of the Acoustical Society of America 131 (2012) 3470–3470. [CrossRef] [Google Scholar]
 N. Lopes, T. Hélie, R. Caussé: Control of an artificial mouth playing a trombone and analysis of sound descriptors on experimental data, in Proceedings of the Stockholm Music Acoustics Conference 2013, July 30–August 3, 2013, Stockholm, Sweden. https://hal.archivesouvertes.fr/hal01245388/document, 2013. [Google Scholar]
 J.B. Doc, C. Vergez: Oscillation regimes produced by an alto saxophone: Influence of the control parameters and the bore inharmonicity. The Journal of the Acoustical Society of America 137 (2015) 1756–1765. [CrossRef] [Google Scholar]
 J. Gilbert, S. Maugeais, C. Vergez: From the bifurcation diagrams to the ease of playing of reed musical instruments. A theoretical illustration of the bouassebenade prescription? in Proceedings of the International Symposium on Musical Acoustics, September 13–17, 2019, Detmold, Germany. http://pub.degaakustik.de/ISMA2019/data/articles/000057.pdf, 2019. [Google Scholar]
 V. Fréour, H. Masuda, S. Usa, E. Tominaga, Y. Tohgi, B. Cochelin, C. Vergez: Numerical analysis and comparison of brass instruments by continuation, in Proceedings of the International Symposium on Musical Acoustics, September 13–17, 2019, Detmold, Germany. http://pub.degaakustik.de/ISMA2019/data/articles/000015.pdf, 2019. [Google Scholar]
 R. Tournemenne, J.F. Petiot, B. Talgorn, M. Kokkolaras, J. Gilbert: Brass instruments design using physicsbased sound simulation models and surrogateassisted derivativefree optimization. Journal of Mechanical Design 139 (2017) 041401. [CrossRef] [Google Scholar]
 S.N. Rasband: Chaotic dynamics of nonlinear systems. Courier Dover Publications, 2015. [Google Scholar]
 R. Seydel: Practical bifurcation and stability analysis, Vol. 5. Springer Science & Business Media, 2009. [Google Scholar]
 I.R. Manchester, M.M. Tobenkin, M. Levashov, R. Tedrake: Regions of attraction for hybrid limit cycles of walking robots. IFAC Proceedings Volumes 44 (2011) 5801–5806. [CrossRef] [Google Scholar]
 A. Schwab, M. Wisse: Basin of attraction of the simplest walking model, in Proceedings of the ASME Design Engineering Technical Conference, Vol. 6. 2001, pp. 531–539. [Google Scholar]
 M. Wisse, A.L. Schwab, F.C. van der Helm: Passive dynamic walking model with upper body. Robotica 22 (2004) 681–688. [CrossRef] [Google Scholar]
 M. Wisse, A.L. Schwab, R.Q. van der Linde, F.C. van der Helm: How to keep from falling forward: Elementary swing leg action for passive dynamic walkers. IEEE Transactions on Robotics 21 (2005) 393–401. [CrossRef] [Google Scholar]
 T. Matsumoto: A chaotic attractor from chua’s circuit. IEEE Transactions on Circuits and Systems 31 (1984) 1055–1058. [CrossRef] [Google Scholar]
 G. Pegna, R. Marrocu, R. Tonelli, F. Meloni, G. Santoboni: Experimental definition of the basin of attraction for chua’s circuit. International Journal of Bifurcation and Chaos 10 (2000) 959–970. [CrossRef] [Google Scholar]
 N.V. Stankevich, N.V. Kuznetsov, G.A. Leonov, L.O. Chua: Scenario of the birth of hidden attractors in the chua circuit. International Journal of Bifurcation and Chaos 27 (2017) 1730038. [CrossRef] [Google Scholar]
 T. Idogawa, T. Kobata, K. Komuro, M. Iwaki: Nonlinear vibrations in the air column of a clarinet artificially blown. The Journal of the Acoustical Society of America 93 (1993) 540–551. [CrossRef] [Google Scholar]
 K. Takahashi, H. Kodama, A. Nakajima, T. Tachibana: Numerical study on multistable oscillations of woodwind singlereed instruments. Acta Acustica United with Acustica 95 (2009) 1123–1139. [CrossRef] [Google Scholar]
 S. Terrien, C. Vergez, B. Fabre: Flutelike musical instruments: A toy model investigated through numerical continuation. Journal of Sound and Vibration 332 (2013) 3833–3848. [CrossRef] [Google Scholar]
 S. Brezetskyi, D. Dudkowski, T. Kapitaniak: Rare and hidden attractors in van der polduffing oscillators. The European Physical Journal Special Topics 224 (2015) 1459–1467. [CrossRef] [EDP Sciences] [Google Scholar]
 L. Wang, W. Lu, T. Chen: Multistability and new attraction basins of almostperiodic solutions of delayed neural networks. IEEE Transactions on Neural Networks 20 (2009) 1581–1593. [Google Scholar]
 W.L. Coyle, P. Guillemain, J. Kergomard, J.P. Dalmont: Predicting playing frequencies for clarinets: A comparison between numerical simulations and simplified analytical formulas. The Journal of the Acoustical Society of America 138 (2015) 2770–2781. [Google Scholar]
 E. Petersen, P. Guillemain, J. Kergomard, T. Colinot: The effect of the cutoff frequency on the sound production of a clarinetlike instrument. The Journal of the Acoustical Society of America 145 (2019) 3784–3794. [Google Scholar]
 A. Guilloteau, P. Guillemain, J. Kergomard, M. Jousserand: The effect of the size of the opening on the acoustic power radiated by a reed woodwind instrument. Journal of Sound and Vibration 343 (2015) 166–175. [Google Scholar]
 T. Colinot, L. Guillot, C. Vergez, P. Guillemain, J.B. Doc, B. Cochelin: Influence of the “ghost reed” simplification of the bifurcation diagram of a saxophone model. Acta Acustica United with Acustica 105 (2019) 1291–1294. [Google Scholar]
 A. Hirschberg, J. Kergomard, G. Weinreich: Mechanics of musical instruments, Chapter 6 (J. Kergomard): Elementary considerations on reedinstrument oscillations. CISM Courses and Lectures 355 (1995) 229–290. [Google Scholar]
 A. Muñoz Arancón, B. Gazengel, J.P. Dalmont, E. Conan: Estimation of saxophone reed parameters during playing. The Journal of the Acoustical Society of America 139 (2016) 2754–2765. [Google Scholar]
 V. Chatziioannou, M. van Walstijn: Estimation of clarinet reed parameters by inverse modelling. Acta Acustica United with Acustica 98 (2012) 629–639. [Google Scholar]
 S. Bilbao, A. Torin, V. Chatziioannou: Numerical modeling of collisions in musical instruments. Acta Acustica United with Acustica 101 (2015) 155–173. [Google Scholar]
 J. Backus: Smallvibration theory of the clarinet. The Journal of the Acoustical Society of America 35 (1963) 305–313. [Google Scholar]
 A. Hirschberg, R. Van de Laar, J. MarrouMaurieres, A. Wijnands, H. Dane, S. Kruijswijk, A. Houtsma: A quasistationary model of air flow in the reed channel of singlereed woodwind instruments. Acta Acustica United with Acustica 70 (1990) 146–154. [Google Scholar]
 J.P. Dalmont, J. Gilbert, S. Ollivier: Nonlinear characteristics of singlereed instruments: Quasistatic volume flow and reed opening measurements. The Journal of the Acoustical Society of America 114 (2003) 2253–2262. [Google Scholar]
 J.P. Dalmont, J.C. Le Roux: A new impedance sensor for wind instruments. The Journal of the Acoustical Society of America 123 (2008) 3014. [Google Scholar]
 A. Chaigne, J. Kergomard: Acoustique des instruments de musique (Acoustics of musical instruments). Belin, 2008. [Google Scholar]
 J. Kergomard, P. Guillemain, F. Silva, S. Karkar: Idealized digital models for conical reed instruments, with focus on the internal pressure waveform. The Journal of the Acoustical Society of America 139 (2016) 927–937. [Google Scholar]
 F. Silva, C. Vergez, P. Guillemain, J. Kergomard, V. Debut: Moreesc: A framework for the simulation and analysis of sound production in reed and brass instruments. Acta Acustica United with Acustica 100 (2014) 126–138. [Google Scholar]
 P. Guillemain, J. Kergomard, T. Voinier: Realtime synthesis of clarinetlike instruments using digital impedance models. The Journal of the Acoustical Society of America 118 (2005) 483–494. [Google Scholar]
 T. Colinot: Numerical simulation of woodwind dynamics: Investigating nonlinear sound production behavior in saxophonelike instruments. PhD thesis, AixMarseille Université, 2020. [Google Scholar]
 J.B. Doc, C. Vergez, S. Missoum: A minimal model of a singlereed instrument producing quasiperiodic sounds. Acta Acustica United with Acustica 100 (2014) 543–554. [Google Scholar]
 J. Gilbert, J. Kergomard, E. Ngoya: Calculation of the steadystate oscillations of a clarinet using the harmonic balance technique. The Journal of the Acoustical Society of America 86 (1989) 35–41. [Google Scholar]
 N.M. Krylov, N.N. Bogoliubov: Introduction to nonlinear mechanics. Princeton University Press, 1949. [Google Scholar]
 M. Nakhla, J. Vlach: A piecewise harmonic balance technique for determination of periodic response of nonlinear systems. IEEE Transactions on Circuits and Systems 23 (1976) 85–91. [Google Scholar]
 B. Cochelin, C. Vergez: A high order purely frequencybased harmonic balance formulation for continuation of periodic solutions. Journal of Sound and Vibration 324 (2009) 243–262. [Google Scholar]
 L. Guillot, B. Cochelin, C. Vergez: A Taylor seriesbased continuation method for solutions of dynamical systems. Nonlinear Dynamics 98 (2019) 2827–2845. [Google Scholar]
 B. Bentvelsen, A. Lazarus: Modal and stability analysis of structures in periodic elastic states: Application to the Ziegler column. Nonlinear Dynamics 91 (2018) 1349–1370. [Google Scholar]
 L. Guillot, A. Lazarus, O. Thomas, C. Vergez, B. Cochelin: A purely frequency based floquethill formulation for the efficient stability computation of periodic solutions of ordinary differential systems. Journal of Computational Physics 416 (2020) 109477. [Google Scholar]
 A. Lazarus, O. Thomas: A harmonicbased method for computing the stability of periodic solutions of dynamical systems. Comptes Rendus Mécanique 338 (2010) 510–517. [Google Scholar]
 T. Colinot, P. Guillemain, C. Vergez, J.B. Doc, P. Sanchez: Multiple twostep oscillation regimes produced by the alto saxophone. The Journal of the Acoustical Society of America 147 (2020) 2406–2413. [Google Scholar]
 W. Beyn, A. Champneys, E. Doedel, W. Govarets, U. Kuznetsov, A. Yu, B. Sandstede: Handbook of Dynamical Systems (Vol 2), Chapter Numerical Continuation, and Computation of Normal Forms. Elsevier, 2002. [Google Scholar]
 Y.A. Kuznetsov: Elements of applied bifurcation theory, Vol. 112. Springer Science & Business Media, 2013. [Google Scholar]
 B. Bergeot, A. Almeida, C. Vergez, B. Gazengel: Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure. Nonlinear Dynamics 73 (2013) 521–534. [Google Scholar]
 F. Silva: Émergence des autooscillations dans un instrument de musique à anche simple. PhD thesis, AixMarseille Université, 2009. [Google Scholar]
 P. Guillemain, C. Vergez, D. Ferrand, A. Farcy: An instrumented saxophone mouthpiece and its use to understand how an experienced musician plays. Acta Acustica United with Acustica 96 (2010) 622–634. [Google Scholar]
 A.H. Benade: Fundamentals of musical acoustics. Courier Corporation, 1990. [Google Scholar]
 A.H. Benade, D. Gans: Sound production in wind instruments. Annals of the New York Academy of Sciences 155 (1968) 247–263. [Google Scholar]
 H. Bouasse: Instruments à vent. Impr. Delagrave, 1929. [Google Scholar]
 D.M. Campbell, J. Gilbert, A. Myers: The Science of Brass Instruments. Springer, 2020. [Google Scholar]
 J.M. Chen, J. Smith, J. Wolfe: Saxophone acoustics: Introducing a compendium of impedance and sound spectra. Acoustics Australia 37 (2009) 18–23. [Google Scholar]
Cite this article as: Colinot T. Vergez C. Guillemain P. & Doc JB. 2021. Multistability of saxophone oscillation regimes and its influence on sound production. Acta Acustica, 5, 33.
All Tables
Parameters of the numerical model: musician control parameters γ and ζ, reed parameters q_{r} and ω_{r}, contact parameter K_{c} and parameters inherent to the numerical implementation η, N_{m} and H.
All Figures
Figure 1 Measured input impedance (solid) and modal reconstruction (dashed) for the fingering D♯ of the alto saxophone. 

In the text 
Figure 2 Timedomain synthesized pressure signal. (a) temporal envelope (black) and blowing pressure parameter γ (red). (b) Normalized spectrogram (dB) with regime names indicated (unstable ones between parentheses): equilibrium (Eq.), second register, quasiperiodic and first register. The signal is extracted from the same blowing pressure ramp as in Figure 4, between γ = 0.45 and γ = 0.51 (at first occurrence of oscillation). 

In the text 
Figure 3 Bifurcation diagram obtained with harmonic balance method and numerical continuation: L^{2}norm of the acoustical pressure depending on the blowing pressure parameter γ for the low written D♯ fingering of an alto saxophone. Thick lines: stable solution, thin lines: unstable solutions. Black: equilibrium, green: 1st register regimes, red: 2nd register regimes. Multistability zones are shaded: light yellow where 1st and 2nd register coexist, darker gray for equilibrium and 1st register. Blue circles specify the location of bifurcations. Vertical black lines correspond to those in Figure 6c, and (from left to right) to phase diagrams 7a–7c. 

In the text 
Figure 4 L^{2}norm of the timedomain synthesis signal (dark blue line) for control scenario number one, a ramp of the blowing pressure parameter γ from 0 to 2 and down, overlaid with the bifurcation diagram of Figure 3 (thick stable solutions and thin unstable solutions, of the first register in green and second register in red). Solid dark blue line indicates the first part of the signal (increasing γ) and dashed dark blue line the second (decreasing γ). 

In the text 
Figure 5 Four examples of blowing pressure variations in control scenario number two, as described by Equation (17), for γ_{f} = 0.4 and 0.6 and τ_{g} = ${\tau}_{g}^{\left(1\right)}$ = 10 ms and ${\tau}_{g}^{\left(2\right)}$ = 20 ms. 

In the text 
Figure 6 Classification of the steadystate regimes produced depending on the blowing pressure transient parameters: final value γ_{f} and characteristic time τ_{g}. Multistability zones (deduced from Fig. 3): (a) equilibrium ○ and first register (b) first register and second register (c) all three regime types. Blue triangles indicate quasiperiodic regimes. The horizontal line on graph (b) shows τ_{g}= 1/f_{1}. The vertical lines highlight the γ values of phase diagrams in Figures 7 and 8. 

In the text 
Figure 7 Projection of the attraction basins obtained with control scenario number three. Large dots are initial conditions, small dots are points the synthesis goes through. The dots’ colors indicate the final regime they lead to (black: equilibrium, green: first register, red: second register). Blue lines represent the limit cycles. The small one on the inside is the first register and the one on the outside is the second register. The blowing pressure γ is (a) 1, (b) 1.1, (c) 1.2 (highlighted in Figs. 3 and 6). 

In the text 
Figure 8 Attraction basins (dots) and limit cycles (lines) in a 3D projection of the phase space obtained with control scenario number three for different blowing pressures γ. Black: equilibrium, green: first register, red: second register, blue: quasiperiodic. 

In the text 
Figure 9 Classification of the regimes produced (○ equilibrium, first register, second register, quasiperiodic) depending on control parameters: γ_{f} and ζ. Each rectangle corresponds to a couple (γ_{f}, ζ) and the points inside indicate the regime for each characteristic time τ_{g} (bottom 0.1 ms, middle 3 ms and top 100 ms). One rectangle on graph (b) is annotated as an example. Graphs correspond to two inharmonicities for low written D♯ fingering (a) f_{2} = 2.065f_{1}, the value measured on a real saxophone and (b) f_{2} = 2f_{1}. 

In the text 
Figure 10 Rate of produced regimes (Eq. (19)) for written low B♯ fingering. Green: first register, red: second register, blue: quasiperiodic. Linestyles indicate the characteristic time. Dotted: τ_{g} = 0.1 ms, dashdot: τ_{g} = 3.2 ms, dashed: τ_{g} = 100 ms, solid: averaged rate. An upward triangle marks the maximum first register averaged rate, a downward triangle marks the minimum second register rate. 

In the text 
Figure 11 Rate of first register regime (Eq. (19)) produced for (a) written low D♯ fingering and (b) written low D fingering. Linestyles indicate the characteristic time. Dotted: τ_{g} = 0.1 ms, dashdot: τ_{g} = 3.2 ms, dashed: τ_{g} = 100 ms, solid: averaged rate. An upward triangle marks the maximums of first register rates. 

In the text 
Figure 12 Rate of produced regimes for the lowest fingerings of the alto saxophone (written pitch). Green: first register, red: second register. An upward triangle marks the maximum first register averaged rate, a downward triangle marks the minimum second register rate. Vertical lines mark the measured inharmonicity values. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.