Multistability of saxophone oscillation regimes and its in ﬂ uence on sound production

– The lowest ﬁ ngerings 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. Time-domain 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 ﬁ rst two resonances is varied in order to ﬁ nd 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 ﬁ ngerings show that a slightly positive inharmonicity, close to that measured on a saxophone, leads to ﬁ rst register oscillations for the greatest range of control parameters. A perfect harmonicity (integer ratio between the ﬁ rst two resonances) decreases ﬁ rst 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.


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][2][3][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 self-oscillating 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 limit-cycle (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 time-domain 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 time-domain synthesis in Section 3 (control scenario number 1). Then, in Section 4, a simple test-case 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.

Numerical simulation framework 2.1 Saxophone model
The saxophone model used in this study is comprised of three main elements: a one degree-of-freedom 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 time-domain 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: 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: 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 time-varying control parameter [4] is the dimensionless blowing pressure c: 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 x r , from [32] for the order of magnitude of the contact stiffness K c .

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: where the two parameters of the reed are its angular eigenfrequency x r and its damping coefficient q r , and T. Colinot et al.: Acta Acustica 2021, 5, 33 the contact force is a function of the dimensionless reed opening x + 1 and is taken from [33], 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 g = 10 À3 to avoid non-differentiability at x = À1 (reed closure), The regularization controlled by parameter g 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 g is not necessary for the time-domain synthesis method to function, it is kept for comparison purposes.

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, where f is the dimensionless control parameter accounting for reed opening at rest, w being the effective width of the reed channel and q 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 g as in Equation (6),

The resonator
The input impedance is used to represent the resonator's acoustical response. The dimensionless input impedance Z(x) 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 cross-section of the crook. A cylindrical tube is added by transfer matrix method [38] in post-processing 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, 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, _ p n ðtÞ À s n p n ðtÞ ¼ C n uðtÞ; 8n 2 ½1; N m ; ð12Þ Reðp n ðtÞÞ: ð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 large-scale 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.

Time-domain synthesis
Equations (4), (7) and (12) are discretized using finitedifference approximations for the time-domain derivatives, Table 1. Parameters of the numerical model: musician control parameters c and f, reed parameters q r and x r , contact parameter K c and parameters inherent to the numerical implementation g, N m and H. 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 c ' 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. Quasi-periodic regimes are well-known on saxophone-like instrument models, documented for instance in [2,8,43].

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 i-th variable X i expands to where there are 2H + 1 complex Fourier coefficients X i,h per variable, and x 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, 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 c. More precisely, knowing a solution to the system for a given value of c, the continuation method finds a solution for a slightly higher or lower values of c, 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.cnrs-mrs.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 time-domain 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][50][51]).

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 time-domain synthesis to exhibit how multistability leads to hysteresis. The correspondance between the two methods (the harmonic balance method and time-domain synthesis) is also checked.

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]   [52]. Figure 3 shows the L 2norm of the pressure signal: where T is the period of the signal, and identifies which regime each branch corresponds to. The blowing pressure parameter c spans the interval between 0 and 2. This bifurcation diagram contains branches corresponding to the so-called equilibrium, where no sound is produced, for the lowest and highest c values. The equilibrium at low c corresponds to the musician not blowing hard enough into the instrument to obtain a sound, while at high c equilibrium means that the reed channel remains closed due to the large pressure difference between mouth and mouthpiece. For intermediate c 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 c 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 quasi-periodic 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 period-doubling bifurcations mark the limits of multistability regions between two oscillating regimes. 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 c = 0.3 and the inverse Hopf bifurcation H1 at c = 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 period-doubling bifurcation PD1 (respectively at c = 0.66 and c = 0.79). The Neimark-Sacker bifurcation NS1 mark the destabilization of the second register and the emergence of a quasi-periodic regime (not represented here), sometimes called multiphonics by musicians. The next coexistence zone occurs in the interval between the two period-doubling bifurcations PD1 and PD2 on the second register branch, where a stable double two-step solution [52] emerges. This coexistence zone is not shaded on the figure, as it could represent less of a musical issue, since double two-step 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 period-doubling 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 c where oscillating solutions exist, including arguably crucial c 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 Time-domain synthesis with blowing pressure ramps (control scenario number one) Once multistability zones are identified, time-domain 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 c: in control scenario number one, the parameter c 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 c slope steep enough to limit dynamical bifurcation delays [55]. Figure 4 shows that the synthesized signal starts from c = 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 quasi-periodic regimes before reaching the first register. Once the first register is established, the branch is followed all the way to extinction (around c = 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, c = 2. The blowing pressure c then starts descreasing, and the system stays at equilibrium until Hopf bifurcation H3 is reached, for c ' 1.05. There, the system jumps to the stable second register regime. The second register branch is followed until the period-doubling bifurcation PD2, where the system briefly follows the double two-step branch. The double two-step 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 period-doubling bifurcation point PD1. Then, something rather surprising occurs: the L 2 -norm of the time-synthesized 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 quasi-periodic 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 quasi-periodic regime  becomes unstable (around c = 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 c are very different. Actually, the two paths only coincide in three regions: the lowest and highest c intervals, for which only the equilibrium is stable, and a very small region around c = 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 time-domain 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 time-domain 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 two-step solution between period doubling PD1 and PD2 are found by the time-domain 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.

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 steady-state 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.

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 c 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 c f , during a certain time determined by the parameter s g . The temporal variation of c is given by the expression, which is differentiable infinitely many times. Figure 5 displays four examples of such transients. Other envelopes (sigmoid, sine branch) were tested, only causing small quantitative changes to the results. Figure 6 shows which established regimes appear in time-domain synthesis depending on s g , for final values c 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 time-domain 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 energy-based 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 quasi-periodic regimes, an attribute that can be observed in Figure 2 is used: the fact that the amplitudes of the harmonics of a quasi-periodic 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 quasi-periodic. 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 (s) and the first register ( ). For final values c f between 0.38 and 0.4, the system can converge to both regime depending on the characteristic rising time s c . It is interesting to note that equilibrium is reached for the longest rising times, i.e. the slowest c variation, whereas the oscillating regime is reached for the shortest rising times. This is understandable as a quick c 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 steady-state 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 quasi-periodic regimes ( ). This is the same quasi-periodic regime that appears in time-domain synthesis in Figure 2, which overlays the unstable portion of the second register branch. There is a particular range of characteristic time s g that seems to produce the first register for a larger range of c f . The corresponding values of s g are close to the period of the first register, represented by the horizontal line in 5 (b).
In the last multistability (0.9 c 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 c f region where all three are produced. This can be explained by analyzing the attraction basins (see Sect. 4.2, Fig. 8).

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 c 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 ill-defined. However, after the transient, the blowing pressure stabilizes around its final value c 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 , _ p 2 ), was chosen as a three-dimensional 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, time-domain 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 I 1 2 ½À0:2; 0:2; p I 2 2 ½À2; 2; _ p I 2 2 ½À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 _ 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 c, 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 shows these initial conditions in the plane (p 2 , _ p 2 ), associated with the regimes they lead to, for three values of the blowing pressure parameter c. 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 c 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 c = 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 c = 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 c 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 c = 1.2, where Figure 6c showed more occurences of the equilibrium than for c = 1.1. This is explained by the attraction basin of the equilibrium being largerthere 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 three-dimensional projection of the phase space (p 1 , p 2 , _ p 2 ), at particular values of c 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 c = 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 c f = 0.6). Then, a quasi-periodic attractor appears in Graph 2 (b), for c = 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 spacecompared, 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 quasi-periodic regimes than first register. A similar interpretation can be formulated with regards to graphs 2 (d) and (e), for c = 0.645 and c = 0.72 respectively. For these values of c 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) 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) (c = 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.

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 so-called Bouasse-Benade prescription [58][59][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 Iðs 2 Þ=Iðs 1 Þ 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 quasi-periodic 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. Quasi-periodic regimes are largely considered undesirable in common musical practice, however they are a common issue on the lowest fingerings of the saxophone.

Regime production regions
Expanding on the idea in Figure 6, one can study the produced regime across the two-dimensional parameter space (c f , f), while still varying the characteristic time s g of control scenario number two. Figure 9 shows the classification of obtained regimes for several combinations (c f , 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 f without also varying the properties of the reed x 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 c f and f and three characteristic times s g , for a total of 192 synthesized signals. The range in parameter f 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 f 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.
Focusing on Figure 9a, several features can be described, and recognized from the situations explored in Section 4 with a fixed f. Coexistence regions can be noticed on most of the map, with a given (c f , 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 c = 0.4, the horizontal threshold near f = 0.2 and the extinction threshold around c = 1.1.
Coexistence situations similar to Figure 6b can also be seen on Figure 9a, for example at f = 1.2 and c 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. 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.

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, where regime i can either be first register, second register or quasi-periodic 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 non-oscillating regimes. Figure 10 depicts the rate of each oscillating regime produced depending on inharmonicity for the lowest fingering of the first register, the written B[ which produces the heard note C3 at 131 Hz. On this figure the regimes are counted separately for each characteristic time s 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 quasi-periodic regimes is also displayed on the figure. Note that inharmonicity values around two lead to less quasi-periodic 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 quasi-periodic 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 quasi-periodic regimes. 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 Figure 9. Classification of the regimes produced (s equilibrium, first register, second register, quasi-periodic) depending on control parameters: c f and f. Each rectangle corresponds to a couple (c f , f) and the points inside indicate the regime for each characteristic time s 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 . 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.
For the D] fingering (Fig. 11a), one can see that depending on the chosen increase duration s 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 s 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.

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 Figure 10. Rate of produced regimes (Eq. (19)) for written low B] fingering. Green: first register, red: second register, blue: quasi-periodic. Linestyles indicate the characteristic time. Dotted: s g = 0.1 ms, dash-dot: s g = 3.2 ms, dashed: s 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 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: s g = 0.1 ms, dash-dot: s g = 3.2 ms, dashed: s g = 100 ms, solid: averaged rate. An upward triangle marks the maximums of first register rates. 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 quasi-periodic regime production (see Fig. 10), which are a known issue in low saxophone fingerings.

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 time-domain 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. 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.