Oscillation threshold of woodwind instruments N. Grand * , J. Gilbert † and F. Laloë ‡

Abstract

We give a theoretical study of the nature of the bifurcations occurring at the oscillation threshold of woodwind instruments, or of physical systems obeying similar non-linear equations of motion. We start from the simplest description of the acoustical behavior these instruments, a mathematical model containing two equations only, one of which is linear but includes delays, while the other is non-linear but has no delay, and discuss its predictions concerning the characteristics of the small oscillations. In particular we study the nature of the bifurcation occurring at threshold; if the bifurcation is direct, the amplitude of the oscillations increases progressively when the control parameter exceeds a threshold value; but, if the bifurcation is inverse, very small oscillations are not necessarily stable and the oscillation may jump discontinuously to a finite amplitude. While direct bifurcations correspond better to what naive intuition would expect, the surprising result of our calculations is their occurrence is by no means the general rule. We also discuss the shape (spectral content) of the small oscillations, and show that they do not always become quasisinusoidal in the limit of infinitely small solutions, in contrast with what is often assumed in the literature (Worman rule). Frequency shifts are investigated as well near threshold. More generally, we show how, despite of the simplicity of the equations of motion themselves, the characteristics of the non linearities of the excitator and of those of the resonator combine to produce a variety of possible behaviors which are not necessarily intuitive.

Introduction

The basic understanding of the physical processes which allow permanent oscillations to be sustained in musical instruments is due to Helmholtz [START_REF] Helmholtz | On the sensation of tone[END_REF] and Rayleigh [START_REF] Rayleigh | The theory of sound[END_REF]. These authors described in general terms how energy could be introduced into passive resonators, in synchronization with their natural oscillations (resonances), and how in a permanent regime the acoustical losses could be exactly compensated by these feeding mechanisms. A good illustration of these ideas can be found for instance in the book of Benade [START_REF] Benade | Fundamental of Musical Acoustics[END_REF], who uses simples analogies to explain the mechanism producing synchronization between a linear resonator and a non linear excitator (see in particular § 20.1). Further progress towards a more detailed theory was by Worman in his thesis [START_REF] Worman | Self sustained oscillations in clarinet like systems[END_REF], which contains a detailed set of equations that model the physical system. Worman also made a first attempt to solve the equations in the Fourier domain but, unfortunately, the calculation was not pushed sufficiently far to reach precise conclusions. Probably, at that time, the task was more difficult than now, since the general physical ideas on the generic properties of non linear dynamical systems and their bifurcations were not as spread as they are presently. Various authors have also attempted, for instance, to define threshold conditions for equations of this type, but curiously most of them have not gone beyond considerations on energy conservation (obviously a necessary condition, but not always sufficient) so that they do not draw general conclusions on the nature of the bifurcation. The natural control parameter in the problem is the slope A of the curve giving the acoustical flow of the excitator as a function of pressure. What is often implicitly assumed is that permanent oscillations may exist as soon as A reaches a threshold value, and that the amplitude of the oscillations grows progressively if A is further increased; at the same time, from an initially almost pure sinusoidal shape, the oscillation is supposed to acquire more and more harmonic content. This latter property is sometimes called the "Worman rule", following Benade [START_REF] Benade | Fundamental of Musical Acoustics[END_REF] (see in particular §23.1), while its conditions of validity of the rule are not specified precisely. Of course, these predictions seem physically very natural, and one expects them to be correct, at least in some cases. On the other hand, after all, they should be derived from the equations of motion since more complex scenarios could perfectly take place; indeed, outside of musical acoustics, we know that there are many physical systems which have a complex threshold behavior, for instance jumping discontinuously from rest to a finite oscillation level.

There are nevertheless a few cases where a detailed study of the nature of the oscillations (amplitude, harmonic content) is already given in the literature. Starting from the simplified equations for musical instrument oscillations derived by Mc Intyre et al. in [START_REF] Mcintyre | [END_REF], Maganza et al [6] study a special case where the res-onator has simple properties (a case that we will call "completely degenerate" in section 4) and show that the mathematical solution is then simply given by an iteration, so that easy complete calculations become possible. Actually, not only the oscillation threshold is then known, but also the nature of the oscillations, which turn out to be perfect square wave functions (instead of quasi sinusoidal functions). Moreover, the theory predicts the occurrence of a whole period doubling scenario for the route to the chaos. We therefore have at least one case where, at the price of stringent simplifications concerning the properties of the resonator, a complete series of predictions is obtained, but it turns out that they do not coincide at all with the naive expectations mentioned above! A more general study is therefore necessary to reconcile the two points of view, in order for instance to determine more precisely their respective domain of validity.

The purpose of this article is to reexamine the question of the nature of the oscillations in a more general context, where all particular assumptions on the characteristics of the resonator are released (no degeneracy for instance). The general spirit is not the elaboration of a detailed model including all subtleties and details of a real instrument (dynamics of the reed, local turbulence, etc.). On the contrary we attempt to focus the interest on generic properties which do not depend on the details of the modeling. For this purpose, we use a simplified physical description of a musical instrument that is hoped to contain enough physics to provide reasonably realistic results, while remaining sufficiently concise to allow analytical calculations; the model is strongly inspired of that of [START_REF] Mcintyre | [END_REF], but it may also include various generalizations. The surprising result that we will find is that, in several cases, the nature of the oscillations is actually the opposite of what a naive energy argument would provide: the threshold bifurcation may often be inverse [START_REF] Bergé | L'ordre dans le chaos[END_REF], which implies that the real threshold may be different from that derived from a linear calculation (a local analysis is not sufficient). In the text below, we show how the particular case studied in [6] falls into our general scheme. We also discuss how, if necessary, the dynamics of the reed as well as the displacement flow that it creates may be incorporated in it. A preliminary description of this work can be found in the thesis of one of the authors [START_REF] Grand | Etude du seuil d'oscillation des sytèmes non-linéaires de type instrument à vent[END_REF] as well as in [START_REF] Grand | Influences des non-linéarités sur les caractéristiques des oscillations des instruments à vent[END_REF]; see also the lectures given in Udine by J. Kergomard [START_REF] Kergomard | Mechanics of musical instruments[END_REF].

The organization of this article is the following: in section 2 we give the equations of the model that are used throughout this article, first in the time domain, then in the frequency domain where the non linearities couple together different frequencies; from this, in section 3, we study the properties of small oscillations (infinitesimal solutions of the system), the values of the parameters for which they exist, as well as the oscillation frequency shifts; then, in section 4, we treat a few interesting particular cases, especially resonators having two or more impedance peaks perfectly harmonic and with exactly the same height (degenerate resonators), which allows us to make contact with the previous work of [6].

The equations of the model 2.1 Notation

A basic reference on the modeling of the oscillations of musical instruments is reference [START_REF] Mcintyre | [END_REF], which has in fact inspired much subsequent work by various authors. Our "minimal model", actually already used in [6], is taken from this reference and consists of a set of two equations only, corresponding to the scheme shown in figure 1. The first equation characterizes the acoustical properties of the resonator, as seen from its excitation end; the second equation characterizes the properties of the excitator, more precisely the "valve effects" [START_REF] Benade | Fundamental of Musical Acoustics[END_REF] that occur at the location of the reed and the very end of the mouthpiece. The system of equations is:

p(t) = dτ G(τ ) f (t -τ ) f (t) = F [p(t), P s ] (1)
where p(t) is the acoustical pressure inside the resonator under the reed (defined, as usual, as the excess pressure with respect to atmospheric pressure) and f (t) is the acoustical flow at this location; here, the acoustical resonator itself is characterized by an impulse response function G(τ ) but, equivalently, one may characterize it by its Fourier transform, which is the acoustical impedance Z(ω), and rewrite the first equation in the form:

p(ω) = Z(ω)f(ω) (2)
where p(ω) and f (ω) are now the Fourier transform p(t) and f (t). In (1), F is a non linear function which includes the effects of reed stiffness, dynamics of the air flow between the tip of the reed and the mouthpiece, etc. P s is the pressure inside the source (the mouth cavity of the player) which provides energy to the system; we will often omit P s from F , since it is a fixed parameter. For a general discussion of the form of this function, see for instance [START_REF] Benade | Fundamental of Musical Acoustics[END_REF] or [11][12]. Fig. 2 shows a typical shape for the curve representing F ; in the absence of oscillation, the solution of the equations of motion is obtained by intersecting this curve with the straight line of equation p = Z(0)f ; we call p 0 and F 0 the corresponding values of the variables, shown in the figure. Following Worman [START_REF] Worman | Self sustained oscillations in clarinet like systems[END_REF], we assume that F is sufficiently regular to be expanded in a Taylor series around the equilibrium point according to:

F (p) = F 0 + A (p -p 0 ) + B (p -p 0 ) 2 + C (p -p 0 ) 3 + ..... (3)
but, otherwise, we make no assumption on the form of this function. We choose the letter A for the coefficient of the linear term in (3), because this coefficient has the same dimension as an acoustical admittance (inverse of an impedance). This parameter plays the role of a control parameter; if A vanishes, no acoustical flow is provided by the excitator for small oscillations, and therefore no energy; no permanent oscillation can therefore maintain itself. The question is then: what will occur when A increases progressively? The first scenario which comes to mind, already mentioned, is that an infinitely small oscillation will appear at some critical value of A, and grows more and more in size (as well as in harmonic content) for larger and larger values of A; nevertheless, we will see below that other scenarios are perfectly possible.

Stability of the trivial solution

Obviously, the system of equations ( 1) always has as a trivial solution p = p 0 and f = F 0 , which trivially corresponds to no oscillation at all around the equilibrium point. The question is whether or not this solution is stable against small perturbations. What happens if, for instance, the second equation is replaced by:

f (t) = F [p(t), P s ] + δu(t) (4)
where δu(t) is some small perturbation to the acoustical flow (it might for instance arise from microtrubulence of air around the tip of the reed)? As long as the effect of an infinitesimal perturbation remains infinitely small, the trivial solution is stable; but, if this is not the case, the solution becomes unstable.

If the system (1) is linearized, which amounts to replacing F (p) -F 0 by A (p -p 0 ), the scheme of figure (1-b) becomes the classical scheme of a linear feedback loop, so that Bode diagrams and the Nyquist theorem may be applied directly [START_REF] Mees | Dynamics of feedback systems[END_REF] [START_REF] Stephani | Design of feedback control systems[END_REF]. A first step is to plot on a polar diagram the variations, in the complex plane, of the complex impedance Z(ω) as a function of ω, as shown in figure 3; in a second step, it is sufficient to multiply Z(ω) by the slope A to obtain the Nyquist diagram associated with the linearized system. For an acoustical resonator resembling more or less a cylindrical (or conical) resonator, a first resonance is expected to occur, at a frequency ω = ω 1 ; it corresponds to a real value Z(ω 1 ) of the acoustical impedance. For a higher value ω = ω ′ 1 , a first anti-resonance of Z(ω) (a resonance of the admittance A(ω)) will then occur and the impedance becomes real again, but with a much smaller value Z = Z ′ 1 . For a still higher value of the frequency, ω = ω 2 , a second resonance occurs and Z(ω 2 ) = Z 2 , and so on: a series of resonances separated by anti-resonances take place in succession.

The Nyquist theorem states in this case that the static equilibrium point is stable as long as the curve showing the product AZ(ω) does not circle around the point of coordinates 1 and 0 (on the real axis). Therefore, as long as AZ max < 1, where Z max is the largest among all values Z 1 , Z 2 , etc. the static equilibrium point is indeed stable; it becomes unstable when A reaches the value Z max -1 ; when A continues to increase, it remains unstable until A reaches the value Z min -1

at which AZ(ω) > 1 for all values of ω.

We remark in passing that the occurrence of stability for very small, as well as for very large, values of A is not unexpected, since there is a symmetry between small values and large values of this slope. This is because, in the equations, the pressure and the acoustical flow basically play a completely symmetrical role; indeed, the system is invariant under the transformation:

p(t) ⇔ f (t) Z(ω) ⇔ A(ω) F [p(t)] ⇐⇒ F -1 [f (t)] (5)
where F -1 is the inverse function of F ; this requires replacing the slope A by its inverse 1/A and implies, therefore, that very large and very small values of A play a similar role.

We now leave the trivial solution and study non-zero, but infinitesimal, solutions; the calculations will be easier in the Fourier domain of frequencies than in the time domain.

Fourier transform

The system written in (1) contains two equations; the first is linear but with a time delay, the second non-linear but without any delay. As soon as we are interested in a non trivial solution, we can no longer take a linear approximation and a more detailed calculation is required. We now study permanent oscillations of non zero amplitude and assume that the solution of the equations is periodic, with some frequency ω/2π, which allows us to write the Fourier series of the pressure in the form:

p(t) -p 0 = n x n e inωt (6)
Since p(t) is real, we have for any value of n:

x -n = (x n ) * (7)
(the star indicates complex conjugation). From (2), we can directly obtain the Fourier components of the acoustical flow as 1 :

f n = A n x n + δ n,0 F 0 (8)
where the A n are defined for all integer (positive or negative) values of n by:

A n = 1 Z(nω) (9)
The A n 's are nothing but the acoustical admittances at an harmonic "comb of frequencies" that are contained in the periodic solution under study. Incidentally, for the moment, we do not fix the frequency of oscillation, but consider it as a free variable, to be determined later.

Since the first equation of system (1) in already included in [START_REF] Grand | Etude du seuil d'oscillation des sytèmes non-linéaires de type instrument à vent[END_REF], the only equation that we now take into account is the second. Calculating the Fourier transform of the square and the cube of a function, we obtain for every harmonic nω a different equation; for instance, if n = 0, we get:

A 0 -A -Bx 0 -C (x 0 ) 2 x 0 = 2B |x 1 | 2 + |x 2 | 2 + |x 3 | 2 + ... + +3Cx 0 2 |x 1 | 2 + 2 |x 2 | 2 + 2 |x 3 | 2 ... + +3C (x 1 ) 2 x * 2 + 2x 1 x 2 x * 3 ..... + c.c. ( 10
) (in the last line, c.c. is for complex conjugate). In the same way, the second equation is obtained for n = 1 and reads:

(A 1 -A ′ )x 1 = 2B [x 2 x * 1 + x 3 x * 2 + x 4 x * 3 + ...] + +3Cx 1 |x 1 | 2 + 2 |x 2 | 2 + 2 |x 3 | 2 + ... + +3C 2x 0 x * 1 x 2 + 2x 0 x * 2 x 3 + (x 2 ) 2 x * 3 + ... (11)
where the coefficient A ′ is defined 2 as:

A ′ = A + 2Bx 0 + 3C (x 0 ) 2 (12)
The equations for n = 2 is:

(A 2 -A ′ )x 2 = B ′ (x 1 ) 2 + 2B [x 3 x * 1 + x 4 x * 2 + ...] + +3Cx 2 2 |x 1 | 2 + |x 2 | 2 + 2 |x 3 | 2 + ... + +3C 2x 0 x * 1 x 3 + 2x 1 x * 2 x 3 + (x 3 ) 2 x * 4 + ... (13)
where:

B ′ = B + 3Cx 0 (14)
while the equations for n = 3 and n = 4 are:

(A 3 -A ′ )x 3 = 2B [x 1 x 2 + x 3 x 0 + x 4 x * 1 + x 5 x * 2 + ...] + +Cx 3 1 + 3Cx 3 2 |x 1 | 2 + 2 |x 2 | 2 + |x 3 | 2 + 2 |x 4 | 2 ... + +3C (x 2 ) 2 x * 1 + 2x 0 x 1 x 2 + ... ( 15
)
2 This coefficient is nothing but the slope of the non linear function F around an new equilibrium point for the pressure, displaced by the amount x 0 with respect to the initial equilibrium (in the absence of oscillation). The distinction between A and A ′ applies only if acoustical impedance at zero frequency Z 0 does not vanish; if it does, A 0 becomes infinite and x 0 = 0, so that A ′ and A merely coincide.

and:

(A 4 -A ′ )x 4 = B ′ (x 2 ) 2 + 2Bx 1 x 3 + 2B [x 5 x * 1 + x 6 x * 2 + ...] + +3Cx 4 2 |x 1 | 2 + 2 |x 2 | 2 + 2 |x 3 | 2 + |x 4 | 2 + ... + +3C (x 1 ) 2 x 2 + (x 2 ) 2 x 0 + 2x 0 x 1 x 3 + 2x * 1 x 2 x 3 + ... (16)
Of course, the series of equations continues beyond n = 4: the non linear equation, initially in the time domain, is equivalent to an infinite system of equations in the frequency domain. Notice that, while A, A ′ , B and C are real coefficients (Taylor expansion of a real function), the A n 's are in general complex. We have not written the equations corresponding to the negative frequencies, since they are merely the complex conjugate of those for positive frequencies, and therefore do not contain any additional information.

Small oscillations

When no oscillation is taking place, all x n 's vanish and all equations are trivially satisfied. We now study the situation where an oscillation of very small amplitude is sustained; what we have to do, mathematically, is to determine the properties of the non trivial solutions of the infinite system of equations, assuming that they are infinitely small.

A first calculation

Among all coefficients x n , some at least must have non zero values for a non trivial solution to be obtained. Let us assume that x 1 is not vanishing and consider this coefficient as a first order quantity; we then assume that all other Fourier coefficients are at least of the same order, and examine which order they actually are (we come back to this assumption in section 3.2.2). In equation [START_REF] Nederveen | Acoustical aspects of woodwind instruments[END_REF], we notice that all terms on the right hand side are at least second order, if not more; on the left hand side, for the orders to match, it is therefore necessary that (A 1 -A ′ ) should be infinitely small as well; this, in turn implies that (A 1 -A), which differs from the preceding quantity only by higher order corrections, is also infinitely small. Consequently, if we choose the parameter A as a control parameter, the critical value is given by:

A = A 1 (17)
This, in passing, shows that at threshold, A 1 is necessarily real, and provides a condition for the possible values of the frequency; we nevertheless postpone this discussion for the moment. If, in addition, we assume that A 1 is different from all the other A n 's (non degeneracy of the resonator), it is easy to see that all other x n 's are necessarily of higher order; for instance, the fact that (A 0 -A) does not vanish at threshold inside the left hand side of [START_REF] Kergomard | Mechanics of musical instruments[END_REF] shows that x 0 is second order in x 1 ; the fact that (A 2 -A) does not vanish either shows that x 2 has the same property; in the same way, the inspection of [START_REF] Manneville | Structures dissipatives, chaos et turbulence[END_REF] shows that x 3 is third order in x 1 , while x 4 fourth order, and so on.

Knowing that, we can solve the system of equations up to third order in x 1 . To second order in this quantity, (13) provides:

x 2 ≃ B (x 1 ) 2 A 2 -A (18)
(to this order, the distinction between A and A ′ in the denominator, or B and B ′ in the numerator, is irrelevant) while, to the same order, (11) provides:

x 0 ≃ 2B |x 1 | 2 A 0 -A (19)
We can now insert these results into [START_REF] Nederveen | Acoustical aspects of woodwind instruments[END_REF] and obtain an equation in x 1 :

A 1 -A -2B 2B |x 1 | 2 A 0 -A x 1 = 2B B (x 1 ) 2 x * 1 A 2 -A + 3C |x 1 | 2 x 1 + 0 (x 1 ) 4 (20)
This is an equation of degree three, which actually reduces to second degree after simplification by x 1 (this is permitted since we are looking for a solution of the system with a non vanishing value of the first Fourier coefficient) and reduces to:

|x 1 | 2 2B 2 A 2 -A + 4B 2 A 0 -A + 3C = A 1 -A (21)
Finally, the calculation provides the result:

|x 1 | 2 = (A 0 -A) (A 1 -A) (A 2 -A) 2B 2 (A 0 + 2A 2 -3A) + 3C (A 0 -A) (A 2 -A) (22)
For the sake of simplicity, in what follows, we may assume that the acoustical impedance of the resonator Z 0 at zero frequency is small3 , and merely replace it by a zero value. This is of course an approximation, but this impedance is indeed small since the continuous impedance of a pipe of a diameter of the order of one centimeter (or more) and of a length of the order of a meter (or less) arises only from weak viscosity effects (Poiseuille flow); if necessary, this approximation could be released without any special difficulty, but just at the price of writing more cumbersome results. When Z 0 = 0, the admittance A 0 becomes infinite, and ( 22) reduces to:

|x 1 | 2 = (A 1 -A) (A 2 -A) 2B 2 + 3C (A 2 -A) (23)
This result, combined with ( 18) and ( 19), gives the characteristics of the small oscillations of the physical system. Obviously, for this oscillation to exist, the right hand side of (23) must be real and positive.

Discussion

Nature of the bifurcation

Equation ( 23) allows us to make predictions on the nature of the bifurcation at threshold, i.e. on what phenomenon takes place when the control parameter A is increased beyond the critical value A = A 1 . For the moment, in order to simplify the physical discussion as much as possible, we only discuss the case where A 2 is also real (we have already seen that, at the frequency of oscillation near threshold, A 1 has to be real), in other words the case where the two resonance frequencies are perfectly harmonic ; this condition is of course not necessarily fulfilled for any resonator and this motivates, in section 3.3, a more general treatment. With this simplification, it is simple to write the criterion for a "direct bifurcation" [START_REF] Bergé | L'ordre dans le chaos[END_REF] to occur (a direct bifurcation takes place when small oscillations exist for values of the control parameter for which the trivial solution is instable; if, on the opposite, the trivial solution is stable, the bifurcation is called "inverse"). If we assume that A 1 corresponds to the highest impedance peak, that is if the small oscillations that we are studying correspond to the lowest threshold4 (in terms of A values), we know from the study of section 2.2 that the trivial solution |x 1 | = 0 is stable when A < A 1 , unstable for higher values of A. Under these conditions, the positivity of the right hand side of (23) when (A -A 1 ) > 0 provides the inequality:

(A 2 -A 1 ) 2B 2 + 3C (A 2 -A 1 ) < 0 (24)
In the plane of the variables (A 2 -A 1 ) and C, this corresponds to two regions limited by two boundaries, a hyperbola and a straight line (the horizontal axis), as shown in figure 4 where the letter D is used for "direct", the letter I for "inverse". The two other regions of the plane correspond to an "inverse bifurcation", a case where the oscillation vanishes instead of appearing when A increases through the value A = A 1 . Figure 4 is somewhat analogous to a phase diagram and, in what follows, we will sometimes refer to it in these terms. Condition (24) may be also written as a condition on coefficient C:

C < - 2B 2 3 (A 2 -A 1 ) (25)
In figure 5, we represent schematically the two possible bifurcations in a different way by plotting the amplitude of the oscillation (characterized, for instance, by |x 1 |) as a function of the control parameter A. A direct bifurcation (region D in figure 4) corresponds to an amplitude of oscillations which grows progressively when the control parameter A increases from the threshold value A = A 1 . In this case, the general properties of the Hopf bifurcation [START_REF] Bergé | L'ordre dans le chaos[END_REF][15] [START_REF] Allwright | Harmonic balance and the Hopf bifurcation[END_REF] imply that the small oscillations are stable. On the other hand, an inverse bifurcation (region I in figure 4) may correspond to a more complicated behavior: the small oscillations, obtained in this case only when the control parameter is smaller than the critical value, that is when the trivial solution is stable, are now unstable [START_REF] Manneville | Structures dissipatives, chaos et turbulence[END_REF][16]; this does not mean that stable oscillations are never possible, but, if they exist, they necessarily occur for larger values of the oscillation amplitude, as symbolized by the dashed lines of figure 5 (their existence depends on higher order terms in the non linear function [START_REF] Manneville | Structures dissipatives, chaos et turbulence[END_REF]). If this is the case, the system now has a threshold value A = A t which is actually smaller than A 1 ; at this value, the oscillation amplitude may jump discontinuously towards a finite level. Clearly, this kind of behavior can not be studied within the perturbative techniques that we have used in this article.

Naively, it would seem natural to expect direct bifurcations to occur in most cases, but figure 4 shows that this is far from being the rule. The interesting part of the figure is actually mostly the upper half, since the lowest part is useful only if we are studying oscillations which do not correspond to the "lowest threshold", as noted above. When A 2 -A 1 is positive, figure 4 shows that, if C = 0 (the acoustical flow is a purely quadratic function of the pressure difference), the bifurcation is always inverse. Direct bifurcation can take place only when the cubic coefficient C has a negative value; the larger the quadratic coefficient B, the smaller the range of direct bifurcations. In other words we see that only some specific class of non linear excitators will give rise to a direct bifurcation, a somewhat unexpected result 5 .

At this point, some comparison between our preliminary results and the literature may be useful. An early reference on the subject is the article by Benade and Gans [17], who start from equations which are similar to ours, but do not derive precise mathematical conclusion on the existence, or non existence, of small oscillations. As already noted in the introduction, the work of Worman [START_REF] Worman | Self sustained oscillations in clarinet like systems[END_REF] is actually a precursor to the present study since it contains accurately all equations in the Fourier domain that are our starting point, but unfortunately without any precise study of the threshold (no inverse bifurcation). In more recent literature, one can nevertheless find equations which give qualitative predictions on oscillation threshold, and which are therefore comparable to [START_REF] Bouasse | Instruments à vent[END_REF] and (23). For instance equation (14.27) of the classical book by Fletcher and Rossing [START_REF] Fletcher | The physics of musical instruments[END_REF] has exactly the same structure but, curiously, does not include factors 2 and 3 in the de-nominator. More importantly, the coefficients which appear in this equation are not the Taylor expansion coefficients of a non linear flow function, but defined in (14.20) as coefficients "which are not just numbers but rather operators involving phase shifts that are frequency dependent ", a comment which certainly does not apply to our coefficients. So, no direct comparison of the two analyses is possible; after studying it, the present authors admit that they do not fully understand this section of reference [START_REF] Fletcher | The physics of musical instruments[END_REF]. In an article on the functioning of brass instruments [START_REF] Elliott | Regeneration in brass instruments[END_REF], Elliott and Bowsher obtain a "regeneration condition" which has, as (24), the form of an inequality, but since they do not give a general form to their non linear function, they just obtain one point in the phase diagram of figure 3 and, therefore, do not provide a complete discussion of its features. The adaptation of our calculations to a more specific (and therefore realistic) model of woodwind instruments has been made by Kergomard [10]; in this model, the control parameter is the source pressure P s instead of the slope A, but all coefficients A, B and C are in fact a function of P s . A general mathematical study of the nature of bifurcations in a system of equations similar to (1) can be found in [START_REF] Allwright | Harmonic balance and the Hopf bifurcation[END_REF]; this author excludes as special cases the "degenerate resonators" that we consider below, but treats multidimensional cases that we have not investigated.

To summarize the results of this preliminary study: assuming that the resonator is non-degenerate (non coincidence of the values of impedances at the various harmonics), we find results that agree with the predictions of the "Worman rule" [START_REF] Benade | Fundamental of Musical Acoustics[END_REF]; see equation [START_REF] Fletcher | The physics of musical instruments[END_REF] for instance, which shows explicitly the property in question, while a generalization for higher harmonics is elementary:

x 3 ≃ C (x 1 ) 3 A 3 -A x 4 ≃ B (x 1 ) 4 (A 4 -A) (A 2 -A) B 2 A 2 -A + C 2A 2 + 3A 3 -5A A 3 -A x 5 ≃ etc.. . (26)
while

|x 1 | itself if proportional to |A -A 1 |.
But, to be valid, these results require very specific conditions from the excitator, related to the existence of a cubic correction to the quadratic non-linearity. Somewhat paradoxically, we find that an inverse bifurcation scenario, with a discontinuous jump in the amplitude of oscillation, is exactly as generic as (if not more than) that of a direct bifurcation. These predictions may be related to the results of Idogawa and coll. [START_REF] Idogawa | Non-linear vibrations in the air column of a clarinet artificially blown[END_REF] who find experimentally that artificial blowing of a clarinet leads to sudden transitions between different regimes, but this is only a conjecture and we have not attempted to make a precise connection.

Other types of oscillations

At the beginning of section 3.1 we assumed that, among all x n 's, the coefficient of lowest order is x 1 , all other coefficients being a least of the same order; from this assumption, we derived that they are actually higher order. But we have not yet examined what happens if the non-vanishing coefficient of lowest order is not x 1 , but x 2 for instance; we now study this possibility. From x 2 e 2iωt , the quadratic term in B generates frequencies 4ω and 0; the cubic term in C generates in addition 6ω, 2ω and -2ω, and these terms will act as "source terms" in the infinite system of equations ( 10), ( 11), [START_REF] Manneville | Structures dissipatives, chaos et turbulence[END_REF], etc. As a consequence, we will now find as lowest order expansions:

x 0 ≃ 2B |x 2 | 2 A 0 -A (27)
and:

x 4 ≃ B (x 2 ) 2 A 4 -A (28)
etc. All these harmonics are even. To higher orders, they will in turn combine inside the non linear function to generate higher harmonics, but clearly all of them will still remain even: no "source term" will ever appear for x 1 , x 3 , etc. What we merely find is that the initial "comb of frequencies" is too tight, since one out of two of the frequencies that it contains is actually irrelevant. In other words, we can simply rename x 2 as x 1 , x 4 as x 2 , etc. and ignore the other frequencies, a process which reduces the new analysis to that of the preceding section. No really new phenomenon then occurs; we just have to redefine the value of the fundamental frequency. Finally, the same simplification applies if, instead of x 2 , higher coefficients define the lowest order term. For instance if, instead of x 2 , we had chosen x 3 as a first order coefficient, the same kind of reasoning would show that only harmonics at multiples of three are present, and the same conclusion would be obtained again; redefining x 3 as x 1 reduces the problem to that treated above.

In conclusion of this section, we find that a scenario of "subharmonic buildup" can never take place for a resonator which has no degeneracy, so that the analysis of the preceding section is sufficient. This does not mean that, for degenerate resonators, the situation remains the same; for instance, if A 2 = A 3 , a subharmonic buildup scenario where x 1 is progressively built from x 2 and x 3 is indeed possible, as discussed in § 4.

General case; frequency shifts

Here we treat a more general case, where we do no longer assume that A 1 and A 2 become simultaneously real for some value of the frequency ω/2π. We call ω 0 a value of ω for which A 1 is indeed real, in other words a pole of the imaginary part [A 1 (ω 0 )] I of the function A 1 (ω) :

[A 1 (ω = ω 0 )] I = 0 (29)
We note A 0 1 the (real) value of A 1 (ω) at this frequency:

A 1 (ω 0 ) = [A 1 (ω 0 )] R = A 0 1 ( 30
)
We now assume that the frequency of oscillation remains close to ω 0 , and write the oscillation frequency as:

ω = ω 0 + δω (31)
while the value of the control parameter A, which is close to the threshold value (A 0 1 ) obtained in the last section, is written:

A = A 0 1 + δA (32)
In the same way, we may write the value of A 1 (ω) in the form:

A 1 (ω) = A 0 1 + δA 1 (ω) (33)
while if δω ≪ ω 0 , we have 6 :

δA(ω) = A 0 1 2iQ δω ω 0 ( 34
)
where Q is the dimensionless quality factor of the impedance peak around which the oscillation takes place. As for A 2 (ω), it may a priori have real and imaginary parts, so that we just write it as:

A 2 (ω) = A 0 2 R + i A 0 2 I + 0(δω) (35)
(it turns out that, for our calculation, the correction in 0(δω) is not needed, so that we do not write it explicitly here). Rewriting equation ( 23) with these notations provides, when limiting the calculation to lowest order in δω and δA:

|x 1 | 2 = 2iA 0 1 Q δω ω 0 -δA A 0 21 + iA 0 2 2B 2 + 3C A 0 21 + iA 0 2 (36)
where we have introduced the notation:

A 0 21 = A 0 2 R -A 0 1 ( 37
)
6 More precisely, the validity of (34) requires that δω ≪ ω 0 /Q. and:

A 0 2 = A 0 2 I (38)
Now, writing that the imaginary part of the right hand side of ( 36) is zero provides the equality:

δω ω 0 = δA 2QA 0 1 × 2B 2 A 0 2 2B 2 A 0 21 + 3C (A 0 21 ) 2 + A 0 2 2 (39)
which relates the frequency shift δω to the excess value δA of the control parameter; if C = 0 (purely quadratic non linear function), this relation simplifies into:

δω ω 0 = A 0 2 A 21 × δA 2QA 0 1 ( 40
)
We can also insert (39) into (36) in order to obtain the dependence of the square of the amplitude of oscillation as a function of δA; after a little algebra we find:

|x 1 | 2 = -K 2B 2 A 0 21 + 3C A 0 21 2 + A 0 2 2 δA ( 41
)
where K is a positive numerical coefficient which plays no role in our discussion.

We have therefore obtained several results. First, equation (41) shows that the results of section 3.2 remain basically valid; the bifurcation is direct if:

2B 2 A 0 21 + 3C A 0 21 2 + A 0 2 2 < 0 (42)
The only difference with (24) is that the coefficient of 3C is now increased by the term A 0 2

2

; as a consequence, the "phase diagram" of figure 4 is modified, and the border between the domains of direct and inverse bifurcations is no longer made of two separate lines, but of only one single cubic curve, as shown in figure 6. We also see in this figure that some inharmonicity of the acoustical resonator favors the range of direct bifurcations; but, as in section 3.2, in the case where the (real) impedance at ω is larger than that at 2ω, negative values of the cubic coefficient C are still required for getting a direct bifurcation. The second result it that, as soon as this inharmonicity occurs, a frequency shift appears, which increases linearly with the square of the amplitude of oscillation; the amplitude of the shift is inversely proportional to the quality factor Q of the impedance peak at ω, a natural result since a very narrow first resonance peak is likely to determine the frequency more accurately than a broad resonance; the sign of the shift depends of that of A 21 and of C. The physical origin of the shift is nothing but a "frequency pulling" effect of the second resonance peak onto the first, in analogy with similar effects occurring in other oscillators in physics (lasers for instance); see [START_REF] Dalmont | Some aspects of tuning and clean intonation in reed instruments[END_REF], section 4.2 for an experimental study on single reed instruments.

Generalizations

In this section we show that the "minimal model", which we have used as a starting point of all our calculations, may easily be generalized to include more physical effects than initially. In fact, one may directly apply all conclusions of the preceding sections to more complicated situations, just at the price of simple redefinitions of the parameters of the model. In this article, we only discuss one example: the dynamics of the reed and/or the air flow passing under it, but in ref. [START_REF] Grand | Etude du seuil d'oscillation des sytèmes non-linéaires de type instrument à vent[END_REF] more examples are treated explicitly. At the end of this section, in order to also give an idea of the limits of the model, we also discuss another generalization that can not be incorporated as easily into the equations, but requires more substantial rewriting of the algebra.

Including reed dynamics in the model

In the initial model, we have assumed that the acoustical air flow is a function, with no delay, of the pressure difference across the reed. This means that the reed itself follows instantaneously this pressure difference and that, in turn, the hydrodynamic flow of air under the reed also adapts in the same way. In other words, inertia effects are ignored as well in the reed and in the air flow that it determines. Of course, this is an approximation, and one could for instance try and improve the model by adding one more independent variable, the position y(t) of the reed; y(t) may for instance be defined as the distance between the tip of the reed and that of the mouthpiece, so that at any time it corresponds to the width of the channel that is available to the flow of air. The latter may then be taken equal to:

f (t) = F [y(t), P s ] ( 43
)
where P s is the (constant) pressure inside the mouth of the player. We can now close the system of equations by relating the changes in time of the position y(t) to the pressure p(t); we do this by writing a differential equation:

µ d 2 dt 2 y(t) + γ d dt y(t) + ξy(t) = Sp(t) (44)
which introduces new physics since the position y(t) now follows the pressure variation with some phase shift; µ is the effective mass of the vibrating part of the reed, S its area, γ a coefficient characterizing its damping; ξ depends on the stiffness of the reed. The exact type of differential equation that we choose is irrelevant for the present discussion; here, in order to make things more concrete, we select the familiar case of a second order differential equation, but a differential equation of higher order would be treated exactly in the same way (provided it remains linear). Together with the second line of (1), or with (2), equations ( 43) and (44) form a system of three coupled equations, which a priori may have a more general behavior at threshold than (1). We will see that this is not the case.

To show this, it is sufficient to notice that, since the Fourier version (44) is merely:

p(ω) = C(ω)y(ω) (45)
where y(ω) is the Fourier transform of y(t), and where:

C(ω) = S -1 µω 2 + iγω + ξ (46)
We may therefore eliminate p by condensing ( 2) and ( 44) into a single equation:

f (ω) = C(ω)A(ω)y(ω) (47)
Finally, instead of taking three independent dynamical variables, p, f and y, we may eliminate the pressure and treat the problem with two variables only. In the new version of the elementary model, the flow f still plays the same role, but the role of the acoustical pressure is now played by the position of the reed y(t); under these conditions, the equations which relate the variables are exactly of the same type, with the only change that, according to relation (47), that the acoustical admittance A(ω) should now be replaced by the product C(ω)A(ω). With this substitution, all reasonings of the preceding sections immediately apply; the only change is summarized in the equation7 :

A n ⇒ A ′ n = S -1 µn 2 ω 2 + iγnω + ξ A n (48)
The role of the cubic term in the non linear function, for instance, is therefore the same as discussed above. We nevertheless note that, since the relevant admittances are now defined in (48), an "effective inharmonicity" is introduced even if the acoustical resonator itself has perfectly harmonic resonances; even if all the A n 's become real for the same frequency ω/2π, this is no longer true in general for the A ′ n 's, because reed inertia effects introduce phase shifts. Under these conditions, formula (39) may be used to evaluate frequency shifts arising from the frequency pulling effects arising from the reed. For instance, if the resonance frequency of the reed was chosen close to 2ω, it is clear that pulling effect of the second peak could be strongly reinforced, a direct consequence of (48); a generalization to a higher harmonic than n = 2 is straightforward.

In conclusion of this section, we see that our preceding results can easily be adapted in order to incorporate more general cases. In the next section, we give an opposite example, where a more substantial adaptation is needed.

Effects that require a new calculation

In order to show at least one example of a physical situation that does not reduce almost trivially to the "minimum model", we now consider a generalization of equation ( 43) in the form:

f (t) = F [y(t), p(t) -P s ] ( 49
)
This is a natural generalization since the flow of air under the reed moves under the influence of the pressure difference P s -p(t) so that, if p(t) is not much smaller than P s , relation ( 49) is more appropriate than (43). If we now keep a relation such as (44) between y(t) and p(t), we see that the non linear function becomes a function of variables which include various phase shifts instead of keeping its instantaneous character. In this case, there is no simple redefinition procedure which brings us back to the minimum model, and we really have to calculate the Fourier series of a product of several different Fourier series to extend the analysis. Of course, there is no special difficulty in doing so, but the equations become more cumbersome. This is done in the appendix of ref [START_REF] Grand | Etude du seuil d'oscillation des sytèmes non-linéaires de type instrument à vent[END_REF] and, for brevity, we do not reproduce the results here. We just mention that the results are somewhat similar formally, while more delicate to interpret physically since more and more physical effects interfere at the same time.

Degenerate resonators and special cases

It is clear that the analysis of the equations that we have given so far is not general: we have always assumed that, when the control parameter A becomes equal to one admittance A n for n = n 0 , it remains different from all the other A n 's. In other words, in the left hand side of the Fourier equations written in § 2.3, we have assumed that only one of the brackets (A -A n ) can vanish at one time; consequently, among all threshold conditions for A, one corresponds to the lowest values and therefore determines the position of the threshold. We now study the case where two threshold conditions degenerate into one, that is where two impedances become real and equal for a particular frequency. We first assume that A 2 = A 1 and examine how the analysis should be modified in this case; we will examine later the case where A 3 = A 1 , which is more realistic for a clarinet, as well as the case where A 2 = A 3 , which is interesting because it allows subharmonic generation. Finally, we will go to the other extreme and assume that an infinite number of A n 's are equal (complete degeneracy), which corresponds to the model used by Maganza and coll. [6].

Degeneracy 1-2

If8 A 2 = A 1 , we expect that x 1 and x 2 will be of the same order, and we take them as a first order quantity. Equations ( 10) then provides, to lowest order:

x 0 ≃ 2B |x 1 | 2 + |x 2 | 2 A 0 -A (50)
while relations ( 15) and ( 16) lead to:

x 3 ≃ 2Bx 1 x 2 A 3 -A (51)
and:

x 4 ≃ B (x 2 ) 2 A 4 -A (52)
All these Fourier coefficients are second order. We now insert these results into ( 11) and ( 13), and obtain, when retaining terms up to third order:

A 1 -A -4B 2 |x 1 | 2 + |x 2 | 2 A 0 -A x 1 = 2Bx 2 x * 1 + 4B 2 A 3 -A x 1 |x 2 | 2 + +3C |x 1 | 2 + 2 |x 2 | 2 x 1 (53)
and:

A 1 -A -4B 2 |x 1 | 2 + |x 2 | 2 A 0 -A x 2 = B (x 1 ) 2 + 2B 2 2 |x 1 | 2 A 3 -A + |x 2 | 2 A 4 -A x 2 + +3C 2 |x 1 | 2 + |x 2 | 2 x 2 (54)
We immediately see that these equations are markedly different from equation [START_REF] Idogawa | Non-linear vibrations in the air column of a clarinet artificially blown[END_REF], where the right hand side was entirely of third order; here we have second order as well as third order. It makes thus sense to limit ourselves to second order and to simplify the non-linear system into:

(A 1 -A) x 1 = 2Bx 2 x * 1 (A 1 -A) x 2 = B (x 1 ) 2 (55)
Assuming that x 1 does not vanish (we are not interested in a trivial solution), we can rewrite the first equation of this system as:

2Bx 2 (A 1 -A) = x 1 x * 1 = e 2iϕ (56)
where ϕ is the phase of x 1 . We therefore obtain:

x 2 = A 1 -A 2B e 2iϕ (57)
which, inserted in the second equation of the system, gives in turn:

x 1 = ± A 1 -A B √ 2 e iϕ (58)
(in these equations, the phase factor e iϕ as well as the ± can be cancelled by a change of the time origin).

Our results are significantly different from those of section 3.1. What is still true is that oscillations occur at frequencies close to a zero of the imaginary part of the admittance (or, equivalently, of the impedance); this is because A itself is real, and because (57) and (58) require that A 1 -A be infinitely small for our calculations to be consistent. But, now, there is no condition on the sign of this difference; oscillations are possible whether it is positive or negative. In other words, when two threshold conditions A = A 1 and A = A 2 coalesce into one, they give rise to a completely different bifurcation scheme where the amplitude is now directly proportional to the excess control parameter and where small oscillations are possible independently of the sign (as before, it is natural to expect that the oscillations are stable only when A > A 1 ). What is the explanation for a behavior which is so different from that obtained in the absence of degeneracy? Mathematically, the crucial point is the presence of the two second order terms 2Bx 2 x * 1 and b (x 1 ) 2 in the right hand side of the equations; in other words, the situation is different because the two waves at ω and 2ω can beat together and, directly, bring non-linear corrections to the other.

The main result of this section, therefore, is that degeneracy helps obtaining direct bifurcations 9 . In addition to good "collaboration" between the impedance peaks (good harmonicity), emphasized by Bouasse [START_REF] Bouasse | Instruments à vent[END_REF] and Benade [START_REF] Benade | Fundamental of Musical Acoustics[END_REF], we see that it is also important that the values of the impedances themselves match properly. Can this fact be useful in practice for designing instruments, since direct bifurcations are presumably desirable musically for the instruments (to allow pianissimo playing)? Of course, in reality the mathematical degeneracy condition will never occur with infinite precision, so that no real resonator will ever satisfy A 1 = A 2 exactly. But this does not mean that our preceding discussion does not apply in practice, since all physical phenomena behave continuously as a function of A 1 -A 2 . If A 1 and A 2 are only slightly different, there will be a very small range of the control parameter where the discussion of section 3.1 (no degeneracy) will apply; but, as soon as the amplitude of oscillations increases slightly, a crossover will occur and the system will behave as if A 1 and A 2 were equal, that is as discussed in this section. In other words, as soon as one has left a very small domain of inverse bifurcation, it will resemble a direct bifurcation.

Degeneracy 1-3

In the same spirit, we now study the effects of a degeneracy between the first and the third impedance peaks; we therefore assume that, for some frequency:

A 1 = A 3 (59)
When x 1 and x 3 are taken as first order quantities, the Fourier equations of section 2.3 provide:

x 0 = 2B |x 1 | 2 + |x 3 | 2 A 0 -A (60)
as well as:

x 2 = B (x 1 ) 2 + 2Bx 3 x * 1 A 2 -A (61)
and:

x 4 = 2Bx 1 x 3 A 4 -A (62)
If we insert these equations into ( 11) and ( 15), we get:

A 1 -A -4B 2 |x 1 | 2 + |x 3 | 2 A 0 -A x 1 = 2B 2   |x 1 | 2 + 2 |x 3 | 2 x 1 + 3x 3 (x * 1 ) 2 A 2 -A + + 2x 1 |x 3 | 2 A 4 -A + 3C |x 1 | 2 + 2 |x 3 | 2 x 1
(63) and:

A 1 -A -4B 2 |x 1 | 2 + |x 3 | 2 A 0 -A x 3 = 2B 2 (x 1 ) 2 + 2x 3 x * 1 A 2 -A x 1 + 2 |x 1 | 2 + |x 3 | 2 A 0 -A x 3 + 2 |x 1 | 2 x 3 A 4 -A + C (x 1 ) 3 + 3C 2 |x 1 | 2 + |x 3 | 2 x 3
(64) Here the beating notes between ω and 3ω do not allow to reproduce directly these initial frequencies; we therefore expect the situation to be more similar to that of a non degenerate resonator. To show this in the general case would probably require a numerical solution of equations ( 63) and (64). Here, in order to make a simple analytical calculation possible, we make additional assumptions: we assume that the impedances of the resonator at all frequencies but ω and 3ω are negligible, in other words that A 0 , A 2 as well as A 4 are infinite (alternatively, we could assume that the coefficient B has a negligible value). Under these conditions, the system reduces to:

(A 1 -A) x 1 = 3C |x 1 | 2 + 2 |x 3 | 2 x 1 (A 1 -A) x 3 = 3C 2 |x 1 | 2 + |x 3 | 2 x 3 (65)
Both equations can then be simplified by eliminating a trivial (zero) solution and we obtain:

|x 1 | 2 = |x 3 | 2 = A 1 -A 9C (66)
(curiously, at this stage of approximation, only |x 1 | and |x 3 | are fixed, not their phase). Indeed, what we now find is strongly reminiscent of the results of §3.2: the bifurcation is direct is the sign of the cubic constant C is negative. The next step would be to assume finite values for A 0 , A 2 and A 4 , and treat the terms in B of ( 63) and (64) to first order; we do not give this calculation here.

Degeneracy 2-3; subharmonic generation

Another interesting situation is degeneracy between the second and the third peak:

A 2 = A 3 = real number (67)

We then examine whether the 2-3 degeneracy allows subharmonic generation, while it was shown in § 3.2.2 to be impossible in the absence of degeneracy. We then obtain:

x 1 = 2B x 3 x * 2 A 1 -A (68)
and similar relations for x 0 and x 5 , which are also second order in x 2 and x 3 . We can then write the relations:

A 2 -A -4B 2 |x 2 | 2 + |x 3 | 2 A 0 -A x 2 = 4B 2 1 A 1 -A + 1 A 5 -A |x 3 | 2 x 2 + +3C |x 2 | 2 + 2 |x 3 | 2 x 2 (69)
as well as:

A 2 -A -4B 2 |x 2 | 2 + |x 3 | 2 A 0 -A x 3 = 4B 2 |x 2 | 2 + |x 3 | 2 A 0 -A + |x 2 | 2 A 1 -A + + |x 2 | 2 A 5 -A x 3 + 3C 2 |x 2 | 2 + |x 3 | 2 x 3 (70)
The situation is then similar to that of § 4.2. If we assume that B = 0, the solutions of the equations are:

|x 2 | 2 = |x 3 | 2 = A 2 -A 9C ( 71
)
and a similar situation occurs (direct bifurcation if C is negative). Inserting this result into (68) directly provides:

|x 1 | 2 = 2B (A 2 -A) 9C (A 1 -A) 2 (72)
(here also, within the approximations that have been made, the phase of the Fourier components remain undetermined; a higher order calculation would be necessary to determine them). One could easily proceed in the same way and obtain higher order corrections in B. We therefore see that small oscillations can indeed happen if two peaks are degenerate, and beat non linearly to create a subharmonic frequency.

Complete degeneracy

After studying the effect of partial degeneracy of two peaks, we now come to the extreme case where an infinite number of peaks are degenerate. More precisely, we will assume that all resonances are perfectly harmonic, with equal values of the impedance separately for the even and for the odd resonances This situation corresponds precisely to that studied by Maganza and coll. [6], since these authors assume that the impedance of the acoustical resonator is given by:

Z(ω) = Z ∞ ε + i tan ωL/c 1 + iε tan ωL/c (73)
Here Z ∞ is the impedance of air in an infinite space without boundaries, L the length of the resonator and c the velocity of sound; a small number ε accounts for weak losses in the model, with no attempt of treating them in a realistic way (assuming ε = 0 is a simple way of avoiding divergencies in the formulas). When ω = πc/4L, which is the resonance condition for this resonator, odd harmonics correspond to an infinite value of the tangent, while even harmonics correspond to a zero value, so that (73) provides:

A 1 = A 3 = A 5 = ... = ε/Z ∞ A 0 = A 2 = A 4 = ... = 1/εZ ∞ (74)
These equations imply an exact degeneracy of all the resonant peaks at ω, 3ω, 5ω, etc. As recalled in the introduction, the interesting feature of this mathematical model is that it allows a complete and exact solution, not only for small oscillations, but also for finite oscillations; this includes domains where they undergo written in the preceding sections is that they are valid for oscillations of any amplitude -not only for infinitely small oscillations. If nevertheless we assume the amplitude to be infinitesimal, the first equation of this system provides:

(A -A 0 ) x 0 + Ba 2 = 0 (90)

If we insert the corresponding value of x 0 into the second equation, we get:

(A -A 1 ) + C - 2B 2 A 1 -A 0 a 2 = 0 (91)
which shows that criterion of direct bifurcation is:

C - 2B 2 A 1 -A 0 ≃ C + 2B 2 A 0 < 0 (92)
where A 0 is given in (74). We notice that, for small oscillations, the amplitude a is proportional to |A -A 1 | while the offset x 0 is, according to (90), proportional to the difference (A -A 1 ) itself. Condition (92) corresponds to the prediction of Maganza et al. [6], obtained by a completely different method (graphical iteration); the condition for direct bifurcation depends in general on the dissipation10 in the resonator, contained in ε (through A 0 ); nevertheless, for weak dissipation, A 0 becomes very large, and the condition that we obtain is merely that the cubic coefficient C should be negative.

The study of finite oscillations is possible by a simple calculation. From the second equation of (89) we obtain:

a 2 = - 1 C A -A 1 + 2Bx 0 + 3C (x 0 ) 2 (93)
which we can insert into the first equation to obtain an equation with x 0 only:

8C 2 (x 0 ) 3 + 9BC (x 0 ) 2 + C (2A + A 0 -3A 1 ) + 2B 2 x 0 + B (A -A 1 ) = 0 (94)
The corresponding curve is a cubic; but a root of equation ( 94) is acceptable as a solution only if it gives a positive value to the right hand side of (93). In figure 7 we therefore show on the same graph the cubic function as well as the parabola giving the value of:

-C 2 a 2 = 3C 2 (x 0 ) 2 + 2BCx 0 + C (A -A 1 ) (95)
Solutions of the equations are the intersection of the cubic with the horizontal axis which fall "inside" the parabola, while for A = A 1 we have the limit situation where this intersection point is exactly on it (figure 7-a). What happens when A increases beyond this value? The intersection point moves to the left part of the axis and the parabola moves upwards, with the consequence that the acceptable domain for x 0 is reduced, so that two different cases are possible: either we have the situation of fig. 7-b where the intersection point is still in the "good" region, which means that the bifurcation is direct; or we have the situation of figure 7-c and the bifurcation is inverse. We have already seen that ( 92) is the criterion for the first situation to occur. In figure 8 we show how an acceptable value of x 0 may appear for a value of A T of A which is inferior to A 1 , providing a threshold value corresponding to the appearance of a finite level of oscillation (first order transition, a case somewhat similar to that shown in figure 4). An interesting feature of the fully degenerate case is that the iteration method of [6] also provides predictions on the stability of solutions. For instance, if the slope of the iterated function is less than one (condition A < A 1 ), the attractor of the iteration consists of just one fixed point; therefore, if (92) is not satisfied, the oscillations obtained for A < A 1 are unstable and tend to decay until static equilibrium is obtained. Moreover, when A > A 1 , the same method shows that at some point, the square wave at frequency ω/2π becomes unstable against a period doubling phenomenon, while the solution at the initial frequency is still a mathematical solution, but unstable. This phenomenon is reminiscent of the subharmonic generation studied in section 4.3; see also [START_REF] Kergomard | Mechanics of musical instruments[END_REF].

Conclusion

The analysis that we have performed gives a whole set of predictions concerning the characteristics of small oscillations. A variety of results has been obtained concerning the nature of the bifurcations (direct or inverse), which depends on the characteristics of both the resonator and of the excitator in a way which is not always intuitive. Numerical simulations and real experiments are possible in order to check these results as well as the stability of the solutions. We have indeed done computer simulations, which are described in [START_REF] Grand | Etude du seuil d'oscillation des sytèmes non-linéaires de type instrument à vent[END_REF], as well as a few experiments; for lack of space, here we summarize the experiments in a few sentences only. The general idea was to obtain reproducible permanent oscillations of an acoustical resonator with a physical system which mimics the effects of the embouchure, the reed and the mouthpiece, while being more easily controllable and reproducible than real playing. For this purpose, a fast valve was made with a commercial "shaker" in order to control the input of air from a pressurized box into a cylindrical resonator, and a microphone was used to measure the acoustical pressure inside the resonator near the excitation point; the corresponding electrical signal was filtered, controlled by a non linear electronic device creating the equivalent of F , and amplified for driving the shaker. In this way a feedback loop was created, partly acoustical, partly electrical, and by changing the characteristics of the non linear feedback loop made it possible to check condition (25) for direct bifurcation, as well as the frequency shift formula (39). In this way it was possible to check that the mathematical predictions of this article are indeed relevant to the behavior of real physical oscillators.

Conceptually, probably the major result of this work is the that, in addition to good harmonicity of the resonance peaks, which relates only to the values of the resonance frequencies, the value of the peak impedance itself is also crucial to the regime of weak oscillations. Changing theses peak values may result in changing a direct bifurcation into an inverse one. Moreover, we have seen in section 4.2 that the equality between the height of the first and the second impedance peak may be useful in obtaining stable small oscillations. If the first resonance peak is low, we have also discussed how oscillations may be stabilized by the creation of subharmonics, a case which is related to the functioning of saxophones. We therefore believe that the preceding analysis brings a general frame for a better understanding of the mechanism by which permanent oscillations can occur in woodwind instruments.

FIGURE CAPTIONS

Fig. 1: Figure (a) shows a schematic representation of a woodwind instrument, including a source of flow a constant pressure P s , an excitator E which acts as a non linear "valve" (the mouthpiece and the reed for a reed instrument), and a linear resonator R (the body of the instrument). Figure (b) symbolizes the mathematical model used in this article: the excitator is characterized by a non linear function F which, applied to the difference P s -p, provides the acoustical flow f entering the resonator, while p and f are linearly related through the acoustical impedance Z of the resonator (or its admittance). Fig. 2: A typical non-linear function giving the entering air flow F crossing the valve of fig. 1 as a function of the acoustic pressure p at the output side of the valve (at the entrance of the resonator). In the absence of oscillation, the equilibrium values are p 0 and F 0 , obtained by intersecting the curve with a straight line of slope A 0 (the admittance of the resonator at zero frequency, sometimes assumed infinite in the text so that this static equilibrium point goes to the vertical axis). When p = P s , the flow changes sign. The self sustained oscillations can only occur around function points where the slope A of the curve is positive, that is when energy is fed into the resonator by the valve. In the text, we discuss the threshold condition A = A 1 , where A 1 is the (real) value of the admittance at some frequency. Fig. 3: Polar diagram showing the variations of the complex impedance of the resonator Z(ω) as a function of ω. The first resonance occurs when ω = ω 1 and corresponds to a real value Z(ω 1 ) = Z 1 of the acoustical impedance; then, for ω = ω ′ 1 , a first anti-resonance of Z(ω) occurs (resonance of the admittance A(ω)) and the impedance is again real, but with a much smaller value Z = Z ′ 1 ; for a still higher of the frequency, ω = ω 2 , a second resonance occurs and Z(ω 2 ) = Z 2 , and so on. Multiplying this complex impedance by the slope A creates the Nyquist diagram associated with the linearized system of equations around static equilibrium; as long as AZ 1 < 1, the curve obtained does not circle around the point +1 and this static equilibrium is stable; it becomes unstable beyond this value until A reaches the value at which AZ(ω) > 1 for all values of ω, a second crossover point at which the equilibrium becomes stable again.

Fig. 4: "Phase diagram" showing the regions D where the bifurcation is direct as well as those regions I where it is inverse. The horizontal axis shows the values of the coefficient C of the cubic term in the non-linear function F , the vertical axis the difference A 2 -A 1 between the admittance at 2ω and that at ω. If the bifurcation really corresponds to the "first threshold", that is to the smallest value of the control parameter A for which small oscillations are obtained, then necessarily A 2 > A 1 , so that only the upper part of the diagram is relevant. This figure is only valid if the second resonance peak is perfectly harmonic with the first; the general case is shown in figure 6. (a) direct bifurcation, where small oscillations begins to occur while A increases and crosses the threshold value A 1 .

(b) inverse bifurcation, where small oscillations disappear in the same situation. In this case, they are expected to be unstable, and the stable oscillations to be of finite amplitude. Another threshold, A = A T , occurs for a smaller value of A and corresponds to a first order transition (the system jumps there towards a finite level of oscillation, shown with dashed lines).

Fig. 5 :

 5 figure is only valid if the second resonance peak is perfectly harmonic with the first; the general case is shown in figure6. Fig. 5: Diagram showing the amplitude |x 1 | of the first harmonic of the oscillation as a function of the control parameter A. The dashed lines show extrapolations to oscillations of finite amplitudes that are not calculated in this article. The trivial solution |x 1 | is stable until the value A = A 1 is reached. Two cases are possible:(a) direct bifurcation, where small oscillations begins to occur while A increases and crosses the threshold value A 1 .(b) inverse bifurcation, where small oscillations disappear in the same situation. In this case, they are expected to be unstable, and the stable oscillations to be of finite amplitude. Another threshold, A = A T , occurs for a smaller value of A and corresponds to a first order transition (the system jumps there towards a finite level of oscillation, shown with dashed lines). Fig 6: Phase diagram in the general case where the second resonance peak is not in perfect harmonic relation with the first; the horizontal axis shows the value of C, while the vertical axis gives the difference A 0 21 between the real values of the admittances at 2ω and ω. When the imaginary part A

Fig 6 : 2 , 2 .Fig. 7 : 8 :

 62278 figure is only valid if the second resonance peak is perfectly harmonic with the first; the general case is shown in figure6. Fig. 5: Diagram showing the amplitude |x 1 | of the first harmonic of the oscillation as a function of the control parameter A. The dashed lines show extrapolations to oscillations of finite amplitudes that are not calculated in this article. The trivial solution |x 1 | is stable until the value A = A 1 is reached. Two cases are possible:(a) direct bifurcation, where small oscillations begins to occur while A increases and crosses the threshold value A 1 .(b) inverse bifurcation, where small oscillations disappear in the same situation. In this case, they are expected to be unstable, and the stable oscillations to be of finite amplitude. Another threshold, A = A T , occurs for a smaller value of A and corresponds to a first order transition (the system jumps there towards a finite level of oscillation, shown with dashed lines). Fig 6: Phase diagram in the general case where the second resonance peak is not in perfect harmonic relation with the first; the horizontal axis shows the value of C, while the vertical axis gives the difference A 0 21 between the real values of the admittances at 2ω and ω. When the imaginary part A 0 2 of the admittance at 2ω is non zero, the range of direct bifurcations is increased (dashed regions) while the border with the region of inverse bifurcation is now one single cubic curve; for comparison, the dashed lines show the same border when A 0 2 vanishes. When A 0 21 = A 0 2 , the critical value of C has its minimal value C m = -B 2 /3 A 0 2 . Fig. 7: Graphical solution of equation (94), with the constraint that the right hand side of (93) must be positive; the intersections of the cubic curve with the horizontal axis give solutions which are valid provided that they fall "outside" of the parabola. When A = A 1 , as in figure (a), both curves intersect the horizontal axis at the origin. When A > A 1 , if the bifurcation is direct, one has the situation of figure (b), while if it is inverse one has the situation of figure (c). Fig.8: In the case shown in figure (7-c), a threshold for a value A = A T < A 1 may be obtained by the construction shown in this figure; it corresponds to a jump of the solution to a finite level of oscillation.

The term δ n,0 F 0 arises because of the term -p 0 = -F 0 /A 0 in the left hand side of (6).

More precisely, and as the calculation will show, we assume here that Z 0 is negligible as compared to the (real) resonance values of the impedance which sustain the acoustical oscillations.

If this is not the case, oscillations at twice the frequency would be possible for a smaller value of the control parameter A.

In[6], a similar result was obtained by a completely different, graphical, method; for the highly degenerate resonator considered in this reference, the existence of small oscillations turns out to be also controlled by the sign of C (see also § 4.4). Our calculation here shows that the property can be generalized to non-degenerate resonators.

If, instead of a second order equation as in (44), we write a higher order equation, the correction beween brackets in (48) becomes a polynomial of higher degree than two, but otherwise all conclusions remain the same.

A 1 being real near threshold, A 2 is therefore also real. In other words, we are dealing here with two harmonic resonances which, moreover, have the same height for their resonance peaks.

If, in figure4, we assume that A 2 = A 1 , the point representing the system falls on the horizontal axis, which is the border between direct and inverse regions; it is therefore not easy to decide which type the bifurcation is.

In[6], formula[START_REF] Grand | Influences des non-linéarités sur les caractéristiques des oscillations des instruments à vent[END_REF] gives the ε dependence of the iterated function, so that the criterion of direct bifurcation is indeed ε dependent.

ACKNOWLEDGMENTS

The authors are grateful to Paul Manneville and Jean Kergomard for useful discussions and friendly comments. Useful improvements were suggested by the two referees, who are thanked for this help. The Laboratoire Kastler Brossel of the Département de Physique de l'ENS is a laboratoire associé au CNRS et à l'Université Pierre et Marie Curie.

period doubling and chaos. Here we examine how these degeneracy conditions affect the solution of the equations in the Fourier domain, [START_REF] Kergomard | Mechanics of musical instruments[END_REF], [START_REF] Nederveen | Acoustical aspects of woodwind instruments[END_REF], [START_REF] Mees | Dynamics of feedback systems[END_REF], etc.

As an intermediate function, we introduce a periodic square function G(t), with period 2π/ω, which has value -1 when -(π/ω) < t < 0, and +1 when 0 < t < (π/ω). Its Fourier coefficients are the y n 's given by:

The first of these equations implies, by a calculation which is similar to that leading to [START_REF] Kergomard | Mechanics of musical instruments[END_REF], [START_REF] Nederveen | Acoustical aspects of woodwind instruments[END_REF], [START_REF] Mees | Dynamics of feedback systems[END_REF], etc., the following relations between the coefficients y n :

while the second provides:

Knowing this, we may assume that p(t) varies according to:

where a is the amplitude of the square oscillation while x 0 is the time independent offset around which the oscillation is centered; both are unknown quantities at this stage. The x n 's are then given by:

But it turns out that, with this value of the x n 's, all equations (10), ( 11), [START_REF] Mees | Dynamics of feedback systems[END_REF], etc. are satisfied term by term provided adequate values of x and x 0 are chosen.

For instance, if we insert (80) into [START_REF] Kergomard | Mechanics of musical instruments[END_REF], using the first relation (77) and the fact that all x n 's are zero for even values of n, we find the relation:

In the same way, relation [START_REF] Nederveen | Acoustical aspects of woodwind instruments[END_REF] combined with the second equation ( 78) and the zero value of all even coefficients, leads to:

Now, if we go further, we can check that no additional relation between on x 0 and a is obtained. For instance, inserting (80) into (13) provides, with the help (twice) of the third relation (77):

which is an identity. In the same way, (15) now becomes, when taking into account the fourth relation (78):

which is merely equivalent to (82). To generalize the reasoning to any Fourier component, it is sufficient to note that, because of (76), relation (79) implies that:

so that, by grouping the constant terms and those in G(t) in (3), we obtain:

with:

Because the function G(t) has only odd Fourier components which, according to (23), all correspond to the same admittance of the resonator, we merely have:

so that the equations of motion reduce to:

These results are indeed equivalent to (81) and (82); they provide necessary and sufficient conditions on the two coefficients x 0 and a for (79) to give a solution of the equations of motion. The major difference with all equations that we have


Henri Gouin
email: henri.gouin@univ-cezanne.fr

L M M T Box

ENERGY OF INTERACTION BETWEEN SOLID SURFACES AND LIQUIDS

We consider the wetting transition on a planar surface in contact with a semi-infinite fluid. In the classical approach, the surface is assumed to be solid, and when interaction between solid and fluid is sufficiently short-range, the contribution of the fluid can be represented by a surface free energy with a density of the form Φ(ρ S ), where ρ S is the limiting density of the fluid at the surface.

In the present paper we propose a more precise representation of the surface energy that takes into account not only the value of ρ S but also the contribution from the whole density profile ρ(z) of the fluid, where z is coordinate normal to the surface.

The specific value of the functional of ρ S at the surface is expressed in meanfield approximation through the potentials of intermolecular interaction and some other parameters of the fluid and the solid wall.

An extension to the case of fluid mixtures in contact with a solid surface is proposed.

Introduction

The phenomenon of surface wetting is a subject of many experiments. They have already been used to determine many important properties of the wetting behavior for liquids on low-energy solid surfaces 1 . To gain a theoretical explanation of the phenomenon of wetting, a generalized van der Waals model is often used 2 .

In the 1950's Zisman 3 developed an experimental method of characterizing the free energy of a solid surface by the measurement of the contact angle with respect to the liquid-vapor surface tension and by changing the test liquid.

More recent measurements done by Li and all 4,5 improved the understanding of this problem. While the contact angle and surface tension are macroscopic quantities, they have their origin in molecular interactions.

In the following, we use a mean-field theory to investigate how the surface energy is related to molecular interactions. The approximation of mean field theory is too simple to be quantitatively accurate. However it does provide a qualitative understanding and allows one to calculate explicitely the magnitude of the coefficients in our model.

In 1977, John Cahn 1 gave simple illuminating arguments to describe the interaction between solids and liquids. His model is based on a generalized van der Waals theory of fluids treated as attracting hard spheres [START_REF] Rowlinson | Molecular theory of capillarity[END_REF] . It entailed assigning to the solid surface an energy that was a functional of the liquid density "at the surface"; the particular form Φ(ρ S ) = -γ 1 ρ S + γ 2 2 ρ 2 S of this energy, where ρ S is the fluid density at the solid wall, is now widely known in the literature and is due to Nakanishi and Fisher [START_REF] Nakanishi | [END_REF] . It was thoroughly examined in a review paper by de Gennes 8 . To account for the wetting behavior of liquids on solid surfaces, one needs to know Φ(ρ S ); a major object of the present paper is to obtain explicit formulas for the coefficients γ 1 and γ 2 , expressing them in terms of the parameters in assumed microscopic interaction potentials. Three hypotheses are implicit in Cahn's picture.

i) For the liquid density to be taken to be a smooth function ρ(z) of the distance z from the solid surface, that surface is assumed to be flat on the scale of molecular sizes and the correlation length is assumed to be greater than intermolecular distances (as is the case, for example, when the temperature T is not far from the critical temperature T c ).

ii) The forces between solid and liquid are of short range with respect to intermolecular distances.

iii) The fluid is considered in the framework of a mean-field theory. This means, in particular, that the free energy of the fluid is a classical so-called "gradient square functional". The point of view that the fluid in the interfacial region may be treated as bulk phase with a local free-energy density and an additional contribution arising from the nonuniformity which may be approximated by a gradient expansion truncated at the second order is most likely to be successful and perhaps even quantitatively accurate near the critical point [START_REF] Rowlinson | Molecular theory of capillarity[END_REF] . Some numerical calculations based on Cahn's model and comparison with experiments can be found in the literature 2,9 . The aim of this paper, as was stated, is to obtain the values of the different coefficients of the expression of energy due to the interaction between a solid wall and a liquid bulk. These values are associated with intermolecular potentials of the liquid and the solid wall. They take into account the molecular sizes and finite range of interactions between the molecules of the liquid and the solid. The value of coefficients γ 1 and γ 2 are positive. Consequently the γ 1 term describes an attraction of the liquid by the solid, and the γ 2 term a reduction of the attractive interactions near the surface. In fact, the calculation leads to a more complex functional dependence than the one given by Cahn, Nakanishi and Fisher or de Gennes 1,[START_REF] Nakanishi | [END_REF]8 . The energy of the wall takes into account not only the value of the density at the wall but also the gradient of density normal to the wall. This should be more accurate when the variation of density is strong with respect to molecular sizes. The density of the liquid ν 1 is chosen depending on three coordinates x, y, z. The direction to the wall is associated with z, and near the wall, the density may be chosen to be a function only of z. In fact such an assumption does not simplify the calculations. In our expression, the energy smoothly depends on x, y and allows continuous variations of density along the wall. The connection with Nakanishi-Fisher expression is not affected by this more general hypothesis. We denote by E S the new expression of the wall energy.

The calculation may be extended to more complex cases: for example to the case of a liquid mixture in contact with a solid wall. We propose a general expression of the energy density at the solid surface. We note that the distribution of the concentrations of the components may be influenced by a solid wall effect.

Description of the solid-fluid interaction in mean field approximation

In regions of a fluid where the density ρ in nonuniform, the intermolecular potentials produce a force on a given molecule that may generate surface tension effects [START_REF] Rowlinson | Molecular theory of capillarity[END_REF][START_REF] Israelachvili | Intermolecular and surface forces[END_REF] . Using classical molecular theory, it is possible to obtain a system of pressure and capillary tension that is mathematically equivalent to the stress tensor of a continuum model [START_REF] Rowlinson | Molecular theory of capillarity[END_REF][START_REF] Rocard | [END_REF] . Such a description does not take into account the possibility of interaction of the liquid with the solid wall (when the distribution of density varies strongly in the direction normal to a solid surface). These effects constitute the subject of the present paper.

The effects are illustrated in Figure 1 and are described below. We consider a flat solid wall. The approach is aimed to describe the interaction of a macroscopic surface with a bulk fluid. The so-called "hard sphere diameter" of the molecules is denoted by σ for the fluid and τ for the solid. Then, the minimal distance between solid and liquid molecules is δ = 1 2 (σ + τ ).

In mean-field theory, we represent by ψ(r) the intermolecular potential between two molecules of the fluid and respectively χ(r) the potential between a fluid molecule and a wall molecule at separation r. The density of molecules at a point depends on its coordinates x, y, z, but the masses m 1 and m 2 of each molecule of the liquid and solid are assumed to be fixed. The energy corresponding to the action of all molecules of the liquid and the solid on a given molecule 1 located in 0 (see Figure 1) is

W 0 = i m 2 1 ψ(r i ) + j m 1 m 2 χ(r j ). ( 1
)
Figure 1: Molecular layer between the liquid and the solid surface.

The first summation (over i) is over fluid molecules (except for molecule 1), and the second summation (over j) is over the wall molecules. Molecule 1 is in the fluid and r i , r j are distances from molecule i or j to molecule 1.

Denoting by ν 1 dω 1 the number of molecules of fluid in the volume element dω 1 and ν 2 dω 2 the number of molecules of solid in volume element dω 2 , the potential energy resulting from the action of all molecules in the medium on molecule 1 located in 0 may be expressed in a continuous way:

W 0 = Ω1 m 2 1 ψ(r)ν 1 dω 1 + Ω2 m 1 m 2 χ(r)ν 2 dω 2 (2)
where Ω 1 and Ω 2 are the domains occupied by the liquid and the solid.

Potential energy W 0 will be summed over all the molecules of the liquid. So, the first integral is counted twice in the previous integral over the domain Ω 1 occupied by the liquid. In such a condition, we have to consider only the expression

1 2 Ω1 m 2 1 ψ(r)ν 1 dω 1 + Ω2 m 1 m 2 χ(r)ν 2 dω 2 (3)
for the potential energy W 0 .

Let ν 1 (x, y, z) be an analytic function in each point M (x, y, z) of the liquid.

In the next derivation we replace at the point (0, 0, 0) ν 1 with its Taylor expansion in variables x, y, z limited to the second order. Such an assumption means that the forces between the solid and the liquid are of a short range. This is similar to the expansion given by Rocard [START_REF] Rocard | [END_REF] . This means that ψ(r) is a rapidly decreasing function of r. Then, it is only necessary to know the distribution of molecules at a short distance from molecule 1. This case reflects the reality when the main force potentials decrease as r -6 with the distance. So, the potentials decrease very rapidly from the solid wall. The distribution of density inside the solid is assumed to be uniform. We write this expansion: If we note ρ 1 = m 1 ν 10 and ρ 2 = m 2 ν 2 the densities in the liquid at point (0, 0, 0) and inside the solid wall, we obtain

ν 1 (x,
W 0 = 2π m 1 ρ 1 ∞ σ r 2 ψ(r)dr + π m 1 ρ 1 ∞ h r(h -r)ψ(r)dr + 2π m 1 ρ 2 ∞ h r(r -h)χ(r)dr + π 2 m 1 ∂ρ 1 ∂z ∞ h r(h 2 -r 2 )ψ(r)dr + π 3 m 1 ∆ρ 1 ∞ σ r 4 ψ(r)dr + π 12 m 1 ∆ T ρ 1 ∞ h (-2r 4 + 3hr 3 -rh 3 )ψ(r)dr + π 6 m 1 ∂ 2 ρ 1 ∂z 2 ∞ h r(h 3 -r 3 )ψ(r)dr (5)
In Eq. ( 5), h represents the distance between the molecule 1 and the solid wall, ∆ is the Laplace operator and ∆ T is the Laplace operator tangential to the wall. Details of the calculation are given in Appendix 1.

Notice that the last two terms may also be written in the form

π 12 m 1 ∆ρ 1 ∞ h (-2r 4 + 3hr 3 -rh 3 )ψ(r)dr + π 4 m 1 ∂ 2 ρ 1 ∂z 2 ∞ h rh(h 2 -r 2 )ψ(r)dr (6)
The energy density per unit volume at point 0 is

E o = ν 10 W 0 .
In fact, this result is independant of the reference point. If we denote E the energy per unit volume at any point M in the liquid, we obtain

E = 2π ρ 2 1 ∞ σ r 2 ψ(r)dr + π ρ 2 1 ∞ h r(h -r)ψ(r)dr + 2π ρ 1 ρ 2 ∞ h r(r -h)χ(r)dr + π 2 ρ 1 ∂ρ 1 ∂z ∞ h r(h 2 -r 2 )ψ(r)dr + π 3 ρ 1 ∆ρ 1 ∞ σ r 4 ψ(r)dr + π 12 ρ 1 ∆ T ρ 1 ∞ h (-2r 4 + 3hr 3 -rh 3 )ψ(r)dr + π 6 ρ 1 ∂ 2 ρ 1 ∂z 2 ∞ h r(h 3 -r 3 )ψ(r)dr (7)
Expression ( 7) yields

E = Kρ 2 1 + α(h)ρ 1 ρ 2 + β(h)ρ 2 1 + γ(h)ρ 1 ∂ρ 1 ∂h + Kb 2 ρ 1 ∆ρ 1 + B ( 8
)
with

K = 2π ∞ σ r 2 ψ(r)dr , Kb 2 = π 3 ∞ σ r 4 ψ(r)dr α(h) = 2π ∞ h r(r -h)χ(r)dr , β(h) = π ∞ h r(h -r)ψ(r)dr γ(h) = - π 2 ∞ h r(h 2 -r 2 )ψ(r)dr B = π 12 ρ 1 ∆ T ρ 1 ∞ h (-2r 4 +3hr 3 -rh 3 )ψ(r)dr+ π 6 ρ 1 ∂ 2 ρ 1 ∂z 2 ∞ h r(h 3 -r 3 )ψ(r)dr
The constant b appearing in the ρ 1 ∆ρ 1 term denotes the covolume of the liquid as in the van der Waals equation [START_REF] Rocard | [END_REF][START_REF] Gouin | [END_REF] .

The corresponding energy of the fluid is

W = Ω1 E dω 1 ( 9
)
where E is given by Eq. ( 8) .

To take into account the kinetic effects like in Rocard [START_REF] Rocard | [END_REF] , p. 392, we must add the terms corresponding to kinetic pressure to the value of W . The first term Kρ 2 1 yields the internal pressure. Now, we calculate the different values of coefficients in Eq. 8 in a special case of London forces.

Calculations of the energy of interaction in the case of London forces

It is now possible to calculate the value of W for the particular interaction potentials. For example, one can take (see Appendix 2)

ψ(r) = - k r n-1 and χ(r) = - µ r n-1 (10)
In the case of London forces one has 10 n=7.

In fact this rough approximation is valid only at short range from the wall (hypothesis ii). Following the calculations in Appendix 2, we obtain that the two integrals in the B term are of the order of 1 h ; moreover,

α(h) = - µπ 6h 3 , β(h) = kπ 12h 3 , γ(h) = - kπ 8h 2 (11)
Integrals associated with Eq. ( 9) are taken over an interval [δ, d] , where d is the range of molecular forces in the liquid and δ is the minimal distance between solid and liquid molecules.

For the same reasons as in Eq. ( 4), the expansion of ρ 1 , ∂ρ 1 ∂h , ∂ 2 ρ 1 ∂h 2 at the solid wall is taken in the form ) S . Now we make the important assumption that the variations of ρ 1 take into account several molecular ranges. Hence, δ d can be considered as a small parameter. It implies that the first derivative of ρ 1 with respect to h is on the order of ρ 1 d , the second derivative of ρ 1 is on the order of ρ 1 d 2 . Then, the B term can be removed from the integration of E. Consequently, in the calculation of the surface energy in Eq. ( 4), the second derivative terms may be removed (but not in the calculation of the bulk energy of the liquid associated with K(ρ 2 1 + b 2 ρ 1 ∆ρ 1 )). This result is in agreement with some considerations by de Gennes 13 .

                 ρ 1 = ρ S + h ∂ρ S ∂h + h 2 2 ∂ 2 ρ S ∂h 2 + • • • ∂ρ 1 ∂h = ∂ρ S ∂h + h ∂ 2 ρ S ∂h 2 + • • • ∂ 2 ρ 1 ∂h 2 = ∂ 2 ρ S ∂h 2 + • • • (12) Terms ρ S ,
Keeping terms up to the first order in Eq. ( 8) yields

E = K(ρ 2 1 +b 2 ρ 1 ∆ρ 1 )+α(h)ρ 2 ρ S +h ∂ρ S ∂h +β(h) ρ 2 S +2hρ S ∂ρ S ∂h +γ(h)ρ S ∂ρ S ∂h or E = K(ρ 2 1 +b 2 ρ 1 ∆ρ 1 )+α(h)ρ 2 ρ S +β(h)ρ 2 S +hα(h)ρ 2 ∂ρ S ∂h + 2hβ(h) + γ(h) 2 ∂ρ 2 S ∂h
Then, the energy of the fluid is

W = W 1 + W 2 , with W 1 = Ω1 K(ρ 2 1 + b 2 ρ 1 ∆ρ 1 )dω 1 W 2 = ρ 2 ρ S Ω1 α(h)dhdxdy + ρ 2 S Ω1
β(h)dhdxdy

+ ρ 2 ∂ρ S ∂h Ω1 hα(h)dhdxdy + ∂ρ 2 S ∂h Ω1 2hβ(h) + γ(h) 2 dhdxdy with +∞ δ α(h)dh = - µπ 12δ 2 , +∞ δ β(h)dh = kπ 24δ 2 +∞ δ hα(h)dh = - µπ 6δ , +∞ δ hβ(h)dh = kπ 12δ +∞ δ γ(h) 2 dh = - kπ 6δ (13)
We note that ρ

1 ∆ρ 1 = ρ 1 div( --→ grad ρ 1 ) = div(ρ 1 --→ grad ρ 1 ) -( --→ grad ρ 1 )
2 and the Stokes formula yields

W 1 = Ω1 K(ρ 2 1 -b 2 ( --→ grad ρ 1 ) 2 )dω 1 + Σ Kb 2 ρ S ∂ρ S ∂h ds (14)
Here Σ notes the surface of the solid wall corresponding to the boundary of Ω 1 . As it was said in paragraph 2, to take into account the kinetic effects, we must add to energy W 1 the terms corresponding to kinetic pressure.

The first term Kρ 2 1 yields the internal pressure, and consequently

ε = α(ρ 1 , η) - Kb 2 ρ 1 ( --→ grad ρ 1 ) 2
denotes the internal specific energy. Term ρ 1 α(ρ 1 , η) is the bulk internal volume energy as a function of ρ 1 and of the specific entropy η in the liquid. Let us note that

Σ Kb 2 ρ S ∂ρ S ∂h ds = Σ - kπ 6σ ∂ρ 2 S ∂h ds
Then, a straighforward calculation yields the energy of the liquid in the well known form

W = Ω1 ρ 1 α(ρ 1 , η) -Kb 2 ( --→ grad ρ 1 ) 2 dω 1 + Σ E S ds
but with

E S = -γ 1 ρ S + γ 2 2 ρ 2 S -γ 3 ∂ρ S ∂h -γ 4 ∂ρ 2 S ∂h (15)
Here, γ 4 ∂ρ 2 S ∂h means 2ρ S γ 4 ∂ρ S ∂h , and

γ 1 = µπ 12 δ 2 ρ 2 , γ 2 = kπ 12 δ 2
(16)

γ 3 = µπ 6 δ ρ 2 , γ 4 = kπ 6 1 σ - 1 8 δ
E S is the form of the special energy to be added at the solid surface to obtain the total energy of the liquid. In expression (15), the γ 1 term (favoring large ρ S ) describes an attraction of the liquid by the solid. The γ 2 term represents a reduction of the liquid/liquid attractive interactions near the surface : a liquid molecule lying directly on the solid does not have the same number of neighbors that it would have in the bulk. The terms with the coefficients γ 3 and γ 4 also describe a reduction of the liquid/liquid attractive interactions due to the lack of molecules of the liquid near the wall (in the case where ∂ρ S ∂h is positive).

Expressions ( 15) and ( 16) do generalize the results by Nakanishi and Fisher 7 : Expression (15) contains terms (the first two) similar to the expression given by Nakanishi and Fisher, but we find additional terms associated with γ 3 and γ 4 . Additional terms are the correction of the two first terms. Our main interest is to obtain values of coefficients of the energy of the wall as a function of the properties of molecules and to take into account variations of density at the wall strong enough with respect to molecular sizes.

Extension to the case of liquid mixtures in contact with a solid wall

Here we propose an extension of the above theory to the case of liquid mixtures. This example is for a binary mixture, but there is no reason one could not include more species in the mixture.

The hypotheses are the same as in the case of one-component liquid. The only difference will come from the interaction of molecules of the two liquids.

An adaptation of the previous calculations yields the following results. The potential energy resulting from the action of all molecules in the medium on molecule 1 of liquid 1 located in 0 is

W 10 = Ω1 m 2 1 ψ 1 (r)ν 1 dω 1 + Ω1 m 1 m 2 ψ 3 (r)ν 2 dω 1 + Ω2 m 1 m 3 χ 1 (r)ν 3 dω 2 (17)
This energy is only for one species in the liquid. To determine the whole energy, one must first sum over all molecules of species 1 and then do an analogous summation over the molecules of species 2. Potential energy W 10 will be summed over all molecules of the liquid mixture. In this way the contribution of the liquid-liquid integral, W 10 , will be counted twice in the previous integral over the domain Ω 1 .

We have denoted by m i , the mass of the molecule of fluid i (i ∈ {1, 2} ), ψ 1 and ψ 2 are the potentials of interaction between molecules of fluid 1 with themselves and among molecules of fluid 2, ψ 3 is the potential of molecular interaction between the two fluids, χ i (i ∈ {1, 2}) are the intermolecular potentials of fluids with the solid wall, and ν i (i ∈ {1, 2, 3} ) denote the number of molecules in liquids and solid wall per unit volume.

Then, following the procedure of section 2, we obtain :

E 10 = ν 10 W 10 = 2π ρ 2 1 ∞ σ1 r 2 ψ 1 (r)dr + π ρ 2 1 ∞ h r(h -r) ψ 1 (r)dr + 2πρ 1 ρ 3 ∞ h r(r -h)χ 1 (r)dr + π 2 ρ 1 ∂ρ 1 ∂z ∞ h r(h 2 -r 2 ) ψ 1 (r)dr + π 3 ρ 1 ∆ρ 1 ∞ σ1 r 4 ψ 1 (r)dr + π 12 ρ 1 ∆ T ρ 1 ∞ h (-2r 4 + 3hr 3 -rh 3 )ψ 1 (r)dr + π 6 ρ 1 ∂ 2 ρ 1 ∂z 2 ∞ h r(h 3 -r 3 )ψ 1 (r)dr + 2 πρ 1 ρ 2 ∞ σ2 r 2 ψ 3 (r)dr + πρ 1 ρ 2 ∞ h r(h -r)ψ 3 (r)dr + π 2 ρ 1 ∂ρ 2 ∂z ∞ h r(h 2 -r 2 ) ψ 3 (r)dr + π 3 ρ 1 ∆ρ 2 ∞ σ2 r 4 ψ 3 (r)dr + π 12 ρ 1 ∆ T ρ 2 ∞ h (-2r 4 + 3hr 3 -rh 3 )ψ 3 (r)dr + π 6 ρ 1 ∂ 2 ρ 2 ∂z 2 ∞ h r(h 3 -r 3 )ψ 3 (r)dr
where σ i (i ∈ {1, 2}) is the diameter of molecule of liquid i.

We may also repeat the same calculation for E 20 associated with the second component of the mixture. As for a simple fluid, we may take into account the kinetic effects and the internal energy of a nonhomogeneous mixture as in Fleming et al 14 . Then, we obtain the following additional energy at the solid surface for the liquid mixture in the form

E S = -γ 11 ρ 1S -γ 21 ρ 2S + 1 2 ( γ 12 ρ 2 1S + γ 22 ρ 2 2S + 2γ 32 ρ 1S ρ 2S ) -γ 13 ∂ρ 1S ∂h -γ 23 ∂ρ 2S ∂h - γ 14 ∂ρ 2 1S ∂h + γ 24 ∂ρ 2 2S ∂h + 2 γ 34 ∂ρ 1S ρ 2S ∂h (18)
All the coefficients γ ij can be calculated explicitly after the particular form of interaction potentials was chosen. For example in the case of London forces, the values of coefficients related to the densities of the two fluids at the surface are

γ 11 = µ 1 π 12 δ 2 1 ρ 3 , γ 21 = µ 2 π 12 δ 2 2 ρ 3 γ 12 = k 1 π 12 δ 2 1 , γ 22 = k 2 π 12 δ 2 2 ( 19
)
γ 32 = k 3 π 24 ( 1 δ 2 1 + 1 δ 2
2

)
where ρ 3 is the density of the solid, µ i are the coefficients associated with potentials χ i , i ∈ {1, 2}, k i is associated with potentials ψ i , i ∈ {1, 2, 3} and

δ i = 1 2 (σ i + τ ), i ∈ {1, 2}
are the minimal distances between the solid and molecules of the two species of the mixture.

Such an expression allows one to estimate the influence of a solid wall on each component of a fluid mixture. Depending on the values and signs of different coefficients γ ij , one can estimate the magnitude of the attraction or repulsion effects due to the wall. Placing the fluid mixture in contact with specially chosen solid walls may be an efficient way to separate constituents of a molecular mixture.

Conclusion

The energy of a fluid in contact with a solid wall contains a contribution from the solid which may be represented by a surface density function. For a flat wall we obtain this expression by taking into account not only the density of the liquid at the surface but also its normal derivative. All the numerical values of the coefficients of the surface energy functional are calculated in terms of the parameters of molecules in solid and liquid. This energy characterizes the behavior of the surface in contact with the fluid in the wetting transition. The method may be extended to the case of nonflat solid surfaces which is important in catalysis chemistry.

Appendix 1.

Calculation of the value of E in Eq. ( 8)

To obtain the formula for E given in Eq. ( 8), we have to calculate the two integrals:

1 2 Ω1 m 2 1 ψ(r)ν 1 dω 1 and Ω2 m 1 m 2 χ(r)ν 2 dω 2 a) Calculation of 1 2 Ω1 m 2 1 ψ(r)ν 1 dω 1
This integral is associated with the energy of interaction between molecules in the liquid. Let us denote by S(0, h), the domain occupied by the sphere centered at (0, 0, 0) and with radius h (see Figure 1). We introduce (r, θ, ϕ) the spherical coordinates associated with the center of the sphere. Then,

1 2 Ω1 m 2 1 ψ(r) ν 1 dω 1 = 1 2 S(0,h) m 2 1 ψ(r) ν 1 dω 1 + 1 2 Ω1-S(0,h) m 2 1 ψ(r) ν 1 dω 1
Let us note than for any integers p, q, r, and any boundary Σ of a sphere centered at (0, 0, 0) and radius r,

Σ
x 2p+1 y q z r ds = 0 and which gives relation (7). Now, we change variable as in Figure 2. The origin of the third coordinate z is placed at 0 ′ . Consequently z = 0 ′ 0 with 0 the position of molecule 1.

If we take now the origin of the z-axis at the solid wall and change the orientation in such a way that 00 ′′ = h, we obtain that z + h is a constant for a given molecule. Hence, we get

∂ρ 1 ∂h = - ∂ρ 1 ∂z , ∂ 2 ρ 1 ∂h 2 = ∂ 2 ρ 1 ∂z 2
Then, (we note similarly the Laplace operator in new coordinates and in old coordinates). So, we obtain Eq. ( 8).

E = 2πρ

2 Ω1r (h 2 -π 12 ∂ 2 ν 10 ∂x 2 +∂ 2 ν 10 ∂y 2 m 2 1 ∞ h 2r 4 + 6 ∂ 2 ν 10 ∂z 2 m 2 1 ∞ h r (h 3 +

 2212214613 r)ν 1 sin ϕdϕ r 2 dr dθThe expansion up to the second order of ν 1 (see Eq. (4)), yields 1 r) ν 1 dω 1 = π ν 10 m r 2 ) ψ(r) dr+ 3r 3 hrh 3 ψ(r) dr + π r 3 ) ψ(r) drwhich is the desired result.

Figure 2 : 2 = m 1 m 2 ν 2 sin ϕ dϕ r 2 dr dθ = 2π m 1 m 2 ν 2 ∞r

 2222 Figure 2: Representation of the changed variable associated with the molecular layer at the solid surface.

+ π 12 ρ 1 ∆ T ρ 1 ∞h-2r 4 + 3hr 3 -π 6 ρ 1 ∂ 2 ρ 1 ∂h 2 ∞ h r(h 3 -

 1143123 rh 3 ψ(r) dr + r 3 ) ψ(r) dr

  ∂ 2 ν 10 ∂x 2 +y 2 ∂ 2 ν 10 ∂y 2 +z 2 ∂ 2 ν 10 ∂z 2 + 2xy

			y, z) = ν 10 + x	∂ν 10 ∂x	+ y	∂ν 10 ∂y	+ z	∂ν 10 ∂z	+	(4)
	1 2	x 2 ∂ 2 ν 10 ∂x∂y	+ 2xz	∂ 2 ν 10 ∂x∂z	+ 2yz	∂ 2 ν 10 ∂y∂z	+...
	where ν 10 , at point (0, 0, 0) . ∂ν 10 ∂x	,	∂ 2 ν 10 ∂x 2 , . . . represent the values of ν 1 and its derivatives

Acknowledgements: I have greatly benefited from generous discussions and correspondence with Professor B. Widom and D r A. E. van Giessen. Without their interest and help, this paper would be never published. The support of PRC/GdR CNES-CNRS 1185 is gratefully acknowledged.

 Appendix 2.

Some remarks about the potential associated with cohesive forces

In molecular theory it is proved by the virial method than the so-called "coefficient a" of the van der Waals equation can be obtained from

The term ϕ(r) is the magnitude of the attractive force between the two molecules of the fluid at the distance r, and N is the Avogadro number. The potential ψ 1 (r) =ϕ(r) dr is assumed to be zero at infinity.

Take the forces in the form

and Eq. (20) yields

Let us introduce ψ(r) such that m 2 1 ψ(r) = ψ 1 (r). The term ψ gives the potential per unit mass. Then,

Here, N m = M denotes the molar mass of the fluid. In case n = 7 (London forces), we get:

In the following, we use:

Consequence: Calculation of E in Eq. ( 8) for London forces. Now, using direct integration, we obtain the values of the coefficients given in Eq. (8). Since now

and n = 7

we obtain