Martin Devaud
email: devaud@ccr.jussieu.fr

Thierry Hocquet
email: thierry.hocquet@upmc.fr

Jean-Claude Bacri

Valentin Leroy
email: valeroy77@yahoo.fr

The Minnaert bubble : a new approach

Keywords: numbers: 43.20.+g, 47.55.D-, 47.55.dd

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Since the pioneering paper by Minnaert 1 , a wealth of literature has been devoted to the acoustic bubble problem. To make an idea of the infatuation about this question from the Acoustics community, one may simply consider that the large review of this topic by Leighton 2 already contains, in its 1997 issue, not less than about 1500 references. Roughly, this topic can be taken up from a twofold view-point. (i ) The application view-point first. An air bubble in water is a powerful sound scatterer, which makes it a deciding actor in numerous acoustic situations. From the murmur of the brooks to the noise of the oceans, from the acoustic fog-impeded propagation of SONAR beams to the turning to account of bubbles in the blood stream as contrast agent in medical ultrasonic imaging, few aspects, if any, of the sound propagation in bubbly liquids have not been put forward or used here or there yet. (ii ) The acoustic bubble is a fascinating object from the theoretical view-point too. It raises indeed several fundamental questions, some of them are also encountered in the electrodynamics domain. As a consequence, well-tried methods and concepts of the latter domain can be adapted to the acoustic bubble. It is precisely our approach in this paper, which is deliberately devoted to a theoretical questioning of some basic statements about the bubble. In this connection, our purpose is inter alia to show that the Minnaert bubble behaving like a harmonic oscillator at the isentropic approximation is in fact a highly nontrivial issue. Aiming at proposing an accurate meaning to this statement, we discuss several aspects of the problem as the mass renormalization of the "naked" bubble, the radiative damping and the external forcing of the bubble. In this respect, we make a few remarks about the role played by water compressibility in the single-bubble Minnaert oscillation problem. With this aim, our paper is wilfully progressive, and organized as follows.

In section II, we consider the radial eigenmodes of an infinite extent of (compressible) water surrounding an air bubble. Considering the amplitude of the extrapressure at the air-water boundary, we find that some of these eigenmodes -more precisely an infinite discrete series of them -present resonances. The fundamental resonance corresponds to a wavelength which is much larger than the bubble's radius, and is found to coincide with the Minnaert resonance. We then briefly sketch the modifications that should be introduced to take capillarity and dissipation into account. We finally consider the experimental case of a finite extent of surrounding water, thus involving reverberation effects at the outer boundary. The bubble is then placed at the center of a spherical cavity. We reexamine the radial eigenmodes and the existence of resonances among them. A spectacular lowering of the fundamental eigenmode's frequency is observed if the radii ratio of the water cavity and the air bubble is not large enough, thus entailing an exotic and original "hot chocolate effect" [START_REF] Crawford | [END_REF] .

In section III, we consider the same problem as in section II, but from the point of view of the mere air bubble.

The role of water is then simply to impose a dynamical boundary condition at the bubble's surface. We derive simple expressions for the pressures undergone and radiated by this surface, and we relate them to the resonances found in section II.

In section IV, we focus on the low frequency domain, more precisely the domain in which sound propagation can be neglected within the air cavity. The motion of air inside this cavity is then described, by means of a single dynamical variable, as a uniform inflation/deflation movement. In this "breathing" approximation, we derive a motion equation for the bubble that turns out to be that of a harmonic oscillator. In this respect, we parallel our acoustic problem with its electromagnetic equivalent and show that the pressure forces exerted by water onto the bubble's surface are quite analogous to the Abraham-Lorentz forces undergone by an accelerated electric charge. We are then led to distinguish between the "naked" bubble, with a light mass (due to air's inertia), and the "dressed" bubble, with a heavy renormalized mass (due to water's inertia). We also discuss the corresponding damping terms accounting for the acoustic radiation as well as those -viscous and thermalaccounting for dissipation.

In section V we address the following question: when the Minnaert bubble pulsates, the surrounding water moves accordingly back and forth; which part of this water motion should then be regarded as "sound" ? We propose and discuss an answer to this question and we seize this opportunity to interpret the mass renormalization term introduced in section IV.

II. RADIAL EIGENMODES OF A PULSATING BUBBLE

A. One bubble in an infinite extent of water

Let us consider, in an infinite extent of water, an air bubble with equilibrium radius R and center located at the origin of coordinates. At equilibrium pressure P 0 and temperature T 0 , air and water have mass densities ρ a and ρ w , and adiabatic compressibilities χ a and χ w . Throughout this paper, we shall use exclusively lagrangian dynamical variables: the element of fluid (air or water) which stood at equilibrium at point r (and under pressure P 0 ) is located at time t at position r + u( r, t) and under total pressure P ( r, t) = P 0 + p( r, t). Moreover, we consider here only radial motions: u( r, t) = u(r, t) e r , where e r = r/r, and p( r, t) = p(r, t). Furthermore, by care of simplicity, we shall assume that both air and water undergo transformations small enough for the linear theory to be valid. We shall make two further simplifications, ignoring (i ) any air-water surface tension effect (then considering that extrapressure p(r, t) is continuous at r = R) and (ii ) any thermal or viscous dissipation (then regarding the small air or water transformations as isentropic). Both of these approximations are discussed in subsection II B. As a consequence, we have p(r, t) = -1 χ div u(r, t).

(2.1a)

Then, c = 1/ √ ρχ standing for the sound velocity in air (c a ) or water (c w ), and using Newton's law

ρ ∂ 2 u ∂t 2 = --→ grad p, (2.1b) the d'Alembert's propagation equation ∂ 2 p ∂t 2 = c 2 ∆ p is recovered.
In this section, we shall consider monochromatic solutions of this equation of the form

p(r, t) = ℜ{p(r)e -iωt }, u(r, t) = ℜ{u(r)e -iωt }, (2.2)
where p(r) and u(r) are complex amplitudes. The d'Alembert's equation then reads (∆+q 2 )p(r) = 0, where q = ω/c is the wave number. Owing to the spherical symmetry, the laplacian operator ∆ is reduced to 1 r d 2 dr 2 r, and we get

d 2 dr 2 + q 2 (rp) = 0 p(r) = 1 r (α sin qr + β cos qr). (2.3)
In the air domain (r R), to avoid singularity ar r = 0, coefficient β should be zero, so that (2.3) reduces to r R : p(r) = p a (r) = A r sin q a r.

(2.4a)

On the other hand, in the water domain (r R), the solution (2.3) can be rewritten, for the sake of convenience, in the form

r R : p(r) = p w (r) = W r sin(q w (r -R) + ϕ), (2.4b) with ϕ ∈ [0, π].
From the above equations (2.4a,b), the displacement amplitude u(r) can be derived using Newton's law (u = 1 ρω 2 dp dr ). Then, writing that extrapressure p and displacement u should be continuous at the airwater boundary (r = R), we get the relations

A sin x a = W sin ϕ, (2.5a) A ρ a (sin x a -x a cos x a ) = W ρ w (sin ϕ -x w cos ϕ), (2.5b)
with x a = q a R and x w = q w R. Thus, for a given angular frequency ω, wave numbers q a and q w -and consequently angles x a and x w -are fixed. Then, dividing (2.5b) by (2.5a), we are left with

1 ρ a (1 -x a cot x a ) = 1 ρ w (1 -x w cot ϕ), (2.6)
from which we get

sin ϕ = x w x 2 w + 1 + ρw ρa (x a cot x a -1) 2 , cos ϕ = 1 + ρw ρa (x a cot x a -1) x 2 w + 1 + ρw ρa (x a cot x a -1) 2 .
(2.7)

As a conclusion, for any angular frequency ω, one can find a monochromatic solution satisfying the boundary conditions (2.5): the {air-water} system has a continuum of eigenmodes.

One nevertheless may wonder whether there exists some particular values of ω for which "something happens" at r = R. To see this, let us fix the amplitude W of the standing wave in water, and look for the values of ω that maximize the pressure amplitude A R sin x a (see (2.4a)) at r = R. Allowing for (2.5a), sin ϕ should then be equal to unity: this means that r = R corresponds to an antinode for the standing pressure wave p w (r) displayed in (2.4b). Then cos ϕ = 0 and, according to (2.7),

1 + ρ w ρ a (x a cot x a -1) = 0 tan x a = x a 1 -ρa ρw . (2.8)
The above equation (2.8) has no analytical solution. However, it can be graphically solved as shown in figure 1 where functions tan x a and x a /(1 -ρa ρw ) are plotted versus x a . In this figure, we have set ρ a /ρ w = 1.293 10 -3 , which corresponds to P 0 = 1 atm.

As displayed by the figure, "regular" intersections of the straight line with the tan x a curve are found for x a ≃ 4.49, 7.73, . . . , very close to the 3π 2 , 5π 2 , . . . , (2n+1)π 2 abscissae of the vertical asymptots of tan x a . In addition to these regular intersections, and due to the fact that the 1/(1 -ρa ρw ) slope of the straight line is slightly larger than unity (i.e. the slope of the tangent to the curve tan x a in x a = 0), there is a "singular" intersection with x a ≃ 0.0623 ≪ 1, as displayed by the zoom in figure 1; it corresponds to the fundamental (radial) resonance of the {air-water} system. Using expansion tan x = x + x 3 3 + 2x 5 15 + o(x 7 ), angle x a for this fundamental resonance can be easily derived:

x a 1 + x 2 a 3 + 2x 4 a 15 ≃ x a 1 -ρa ρw x 2 a ≃ 3ρ a ρ w 1 - 1 5
ρ a ρ w , (2.9a)
associated with the angular frequency

ω 0 = c a x a R ≃ 1 - 1 10
ρ a ρ w ω M ≃ ω M , (2.9b)
where

ω M = 3 ρ w χ a R 2
(2.9c)

. In addition, since the slope 1/ `1-ρa ρw ´of the straight line is slightly larger than unity, a singular resonance (the fundamental one) is found for xa ≃ 0.0623 ≪ 1, as displayed by the zoom (in which the difference between the slopes has been deliberately exagerated, for the sake of visibility).

is known as the Minnaert angular frequency.

The above results deserve a few comments. The resonances found when solving equation (2.8) are a priori of the Pérot-Fabry type: they can be understood in terms of constructive interferences of a radial pressure wave bouncing inside the bubble between its center and its surface. If the air-water interface's radius of curvature (i.e. the bubble's radius) R is large compared to the acoustic wavelengths (i.e. if x a = q a R, x w = q w R ≫ 1), the situation corresponds to the plane wave configuration of a 1-D Pérot-Fabry device. In this case, the determination of the resonances is very simple : a pressure wave originating in r = 0 is dephased by q a R = x a when it reaches the bubble's surface. After one reflection on this surface (with no dephasing since water's acoustic impedance is larger than air's one) and a there-and-back center ↔ surface additional travel, this same wave is dephased at r = R by 3x a -π ( -π when bouncing in r = 0). Both paths constructively interfere at r = R provided that 3x a -π = x a + 2nπ, i.e. x a = (2n+1)π 2 . As can be checked from figure 1, this is what is roughly found for the "regular" resonances, corresponding to n = 1, 2, 3, . . . (the higher the value of integer n, the smaller the dis-crepancy). But what about the "singular" resonance, for which x a ≪ 1 ? It is hardly credible that such a small dephasing may yield constructive interferences. In fact it does. This paradox is removed by considering the geometry of the problem: since the wavelengths (in air and in water) are no longer small compared to the air-water interface's radius of curvature, the air-water reflection becomes dispersive, involving a dephasing of the pressure wave which is no longer zero but tends towards π in the x a , x w ≪ 1 limit. In this dephasing originates the spectacular shift of the n = 0 resonance from x a = π 2 (as would be the case in 1-D geometry) to x a ≃ 0.06.

In figure 2, we have displayed sin ϕ = R W p(R) as a function of x a , with the same ratio ρ a /ρ w as in figure 1. Considering the 1 √ 2 -linewidth of the peaks (i.e. the half-height linewidth of sin 2 ϕ) defined as δ = x + a -x - a , where x ± a are solutions of the equations (see (2.7)) In this figure 2, the water compressibility has its actual value (χ w ≃ 10 -10 Pa -1 ). If water is assumed to be incompressible, the positions of the peaks are unchanged (equation (2.8) is independent of χ w ), but their widths are then zero: χ w = 0 c w = ∞ δ = 0, hence figure 2 looks like a comb (strictly speaking, the {air-water} system has then a discrete series of eigenmodes). , there is a singular peak for xa ≃ 0.06 ≪ 1 corresponding to the Minnaert resonance. The compressibility of water is taken into account and the peaks have a finite linewidth, as calculated in the text. If water was assumed to be incompressible, the curve would just be a comb.

1 + ρ w ρ a (x ± a cot x ± a -1) = ∓ x w , (2.10
It is worth emphasizing that, in the case of the singular Minnaert resonance, we have x a = q a R ≪ 1, so that the air pressure in the bubble is quasi homogeneous:

p a (r) = A r sin q a r ≃ Aq a . (2.12a)
Accordingly, the radial displacement u(r) reads

u(r) = u a (r) = 1 ρ a ω 2 0 dp a dr = A ρ a ω 2 0 r 2 [q a r cos q a r -sin q a r] ≃ -χ a Aq a 3 r, (2.12b)
which corresponds to a simple "breathing" (inflation/deflation) of the bubble. This breathing motion will be the master hypothesis of section IV. Observe by the way that equation (2.1a) is satisfied by approximate formulas (2.12a,b).

To end this subsection, it is interesting to list here the conditions that should be satisfied by two media, respectively labelled "a" (inner) and "w" (outer) for a Minnaert-like resonance to be possible. (i ) First, as displayed by condition (2.8), one should have ρ a ≪ ρ w : a drop of water surrounded by air will not give the Minnaert resonance. (ii ) Second, for the Minnaert angular frequency ω M to be telled from zero, the width of the singular peak (see the zoom in figure 2) should be smaller than the corresponding x a , thus implying δ ≪ x a , or equivalently (see (2.9a) and (2.11b)), 3χ w ≪ χ a . Observe by the way that this latter inequality implies that, at the Minnaert frequency, one has

x w = Rω M /c w = 3χ w /χ a ≪ 1.
In other words, condition (ii ) means that the bubble's radius is much smaller than the acoustic wavelength in both media (a and w).

Let us now, in the next subsection, seriatim sketch how the above results are changed when surface tension and dissipation are considered.

B. Capillarity and dissipation effects

Let σ stand for the air-water surface tension, P a (R, t) and P w (R, t) for the (total) pressures at r = R on both sides of the air-water interface. Naturally, we have

P a (R, t) -P w (R, t) = 2σ R + u(R, t) . ( 2

.13)

The equilibrium pressure inside the bubble is then

P ′ 0 = P 0 + 2σ R .
Consequently, the equilibrium mass density ρ ′ a and isentropic compressibility χ ′ a of the air read

ρ ′ a = ρ a 1 + 2σ RP 0 , χ ′ a = χ a 1 + 2σ RP 0 -1 . (2.14)
(Note that the sound velocity c a = 1/ √ ρ a χ a , which depends only on temperature, is unaffected by the capillarity-induced P 0 P ′ 0 change of the bubble's equilibrium inner pressure). With σ = 7 . 10 -2 J.m -2 and under atmospheric pressure P 0 = 10 5 Pa, the ratio 2σ/RP 0 equals unity for R = 1.4 µm: corrections (2.14) can be omitted for bubble's radii R larger than 0.1 mm.

In addition to the equilibrium pressure discontinuity, the air-water surface tension produces at the interface a discontinuity of the extrapressure which reads, at first order in u/R and allowing for (2.13),

p a (R, t) -p w (R, t) = - 2σ R 2 u(R, t). ( 2
+ ρ w ρ a + 2σ R χ ′ a x 2 a (x a cot x a -1) = 0. (2.16)
One observes that the above condition (2.16) substantially departs from (2.8) for small values of x a (i.e. small values of ω). It is thus particularly interesting to examine how the lowest resonance, namely the Minnaert one, is shifted by capillarity. For x a ≪ 1, a straightforward calculation gives, at first order in ratio σ/RP 0 :

x 2 a = 3ρ ′ a ρ w 1 - 2σ 3 χ ′ a R ω 0 = ω M 1 + σ RP 0 3γ -1 3γ . (2.17)
where γ = C P /C V is the specific heats ratio of air. In conclusion, the effects of surface tension are negligible, as well as concerns equilibrium state or dynamics, if condition σ/RP 0 ≪ 1 is satisfied. Henceforth we shall consider bubble radii larger than 0.1 mm and ignore capillarity.

Taking dissipation effects into account is a far more involved affair 4 . Both basic equations (2.1a) and (2.1b) are indeed no longer available and should be completed; moreover, they should be associated with a third equation ruling heat transfer, as (briefly) recalled below.

Let S 0 ( r) d 3 r be the equilibrium entropy of the fluid (air or water) element volume d 3 r located at position r. At time t in the course of motion, the entropy of this same element volume (now located at position r + u( r, t)) reads (S 0 ( r) + s( r, t)) d 3 r. Let T ( r, t) = T 0 + θ( r, t) be its temperature. The entropy field s( r, t) was, by definition, zero in the framework of the isentropic approximation (ii ) in subsection II A; it should now be considered. In order to allow for non isentropicity, straightforward thermodynamical considerations lead to generalize (2.1a) in

p(r, t) = - 1 χ div u(r, t) - s S 0 . (2.18a)
On the other hand, in order to allow for viscosity, equation (2.1b) should be substituted by the Navier-Stokes

equation ρ ∂ 2 u ∂t 2 = --→ grad p + η -→ ∆ ∂ u ∂t + η + ζ 3 -→ grad div ∂ u ∂t , (2.18b)
where η and ζ respectively stand for the so-called first and second viscosity coefficients. Furthermore, the (linearized) heat equation reads

∂s ∂t = Λ T 0 ∆ θ, (2.18c)
where Λ is the fluid's (air or water) thermal conductivity.

Looking for monochromatic solutions of the above set of equations (2.18) of the form (2.2) with in addition

s(r, t) = ℜ{s(r)e -iωt }, θ(r, t) = ℜ{θ(r)e -iωt } (2.19)
is doubtless an interesting issue. It would nevertheless divert us from the purpose of the present paper. Therefore, we postpone this discussion to the end of section IV (devoted to the mere Minnaert bubble) where we summarize Prosperetti's approach 5 and approximations to take dissipation into account, and provisionally ignore the related effects.

C. Finite size effects

In fact, in any laboratory experiment, the bubble is placed in an aquarium, i.e. is immersed in a finite extent of water. It is therefore worth wondering how the above results are affected by such a confinement 6 . It is of course not possible to find an analytic answer to this question for any aquarium shape. It is nevertheless possible, as explained hereafter, in the -not very realisticframework of the radial symmetry.

Let us thus now consider our same air bubble, but surrounded by a finite extent of water bounded by a concentric sphere with radius L (L > R). No extrapressure is applied at the boundary r = L. We are concerned again in the eigenmodes of this {air-finite water} system. The problem is rather simple: since, for 0 < r < L, the equations of motion are unchanged, solutions (2.4) still apply. The extrapressure p w should just vanish at r = L, implying sin[q w (L -R) + ϕ] = 0, with ϕ still given by (2.7). This condition if fulfilled for a discrete series of values ω p of the angular frequency ω satisfying

ω p c w (L -R) + ϕ(ω p ) = pπ, (p = 1, 2, . . . ).
(2.20)

The above ω p can be also regarded as the eigen angular frequencies of the water cavity comprised between the bubble's surface (with a (π + 2ϕ)-dephasing of the pressure wave at r = R, as explained in (3.9a)) and the water's outer boundary (with a π-dephasing of the pressure wave at r = L). Not surprisingly, the ω p are c wdependent. Allowing for (2.7), condition (2.20) can also read (compare with (2.8)) and, at atmospheric pressure, χ a ≃ 7 × 10 -6 Pa -1 , condition (2.22) is satisfied for L R ≫ 2 × 10 5 : an aquarium with the depth of an ocean would be necessary to be sure not to miss any resonance of a 1 mm radius bubble ! Another interesting limit-case is that of "small" aquariums. Let us consider equation (2.20) again and look for the fundamental mode p = 1. We have to determine ω 1 satisfying the condition ω1 cw (L -R) = π -ϕ(ω 1 ). If L -R is small enough, one finds the curious following result. Suppose that ω1 cw (L -R) ≪ π 2 , i.e. that the wavelength (in water) is much larger than 4(L -R). Then ϕ(ω 1 ) ≃ π (and consequently ω 1 is larger than the Minnaert resonance angular frequency displayed in (2.9b), as can be checked from the zoom in figure 2). Equation (2.21) can also be simplified in

1 + ρ w ρ a (x a cot x a -1) = -x w cot(q w (L -R)). (2.21) Since ϕ(ω p ) ∈ [0, π],
c w π L -R ≪ c a R ρ a c a ρ w c w πR L -R ≪ χ w χ a , ( 2
1 + ρ w ρ a (x a cot x a -1) ≃ - R L -R tan x a ≃ x a 1 -ρa ρw L L-R , (2.
ω 1 (L) ≃ ω M (L) = L L -R ω M (> ω M ). ( 2

.24)

The above result will be re-interpreted in section V (see equations (5.6b), (5.7)). Observe that, although it has been derived here assuming

ω1(L-R) cw ≪ π 2 , one has ω M (L) → ω M with L R → ∞.
Nevertheless, it should be kept in mind that ω 1 (L) is but the angular frequency of the fundamental mode of the {bubble-finite water} system. In the very sense precised above in this paper, it is not a resonance: there is no loop of the pressure wave p w (r) at r = R, but rather a quasi-node (since ϕ(ω 1 ) ≃ π); consequently the amplitude rp w (r) of the pressure wave increases like sin ω1 cw (r -R) ≃ ω1 cw (r -R) when one departs from the bubble. Let us recall again that this situation occurs when ω1(L-R) cw

≃ ωM cw L(L -R) < π 2 , i.e. (see (2.9c)) for L R < π 2 χa
3χw ≃ 240 at atmospheric pressure. In this connection, it should be emphasized that experimental setups designed to illustrate the Minnaert resonance will, if miniaturized too much, rather exhibit the {aquarium + bubble} fundamental mode, with the misleading quasi-coincidence displayed in (2.24) of their both angular frequencies. In this respect, it is fascinating how a small amount of air can lower the fundamental frequency of a water cavity. As an illustration, let us take L = 5 cm. Without any bubble, one would have ω 1 = πcw L ≃ 10 5 s -1 . With a R = 1 mm air bubble at the center of the cavity, one has ω 1 ≃ ω M ≃ 2 × 10 4 s -1 , i.e. a lowering of more than two octaves. Although the geometry and the boundary conditions are not the same, a cylindrical mug with height and diameter equal to 10 cm will roughly, when filled up with water, behave like our 5 cm-radius spherical cavity. And a single small 1mmradius air bubble at its center will produce a spectacular non-standard "hot chocolate effect". Let us recall that the "standard" hot chocolate effect, as described by Crawford 3 , relies upon the bubble-induced lowering of the effective Wood velocity 7 in bubbly liquids. A remarkable property of this Wood velocity is that it depends only of the air volume fraction (not of the bubbles' sizes), provided that the sound frequency is much lower than the individual bubbles' Minnaert frequency, which is just not the case here.

We shall henceforth assume that the L R ratio is large enough for the infinite water approximation to be valid. In the next section III, we integrate out the water degrees of freedom and consider the same problem as in this section II, but from the mere bubble's point of view.

III. HOW THE BUBBLE FEELS THE WATER

Since there is no air-water action at a distance, the only outer parameter that rules the dynamics of the air bubble is the water extrapressure p w (R, t) upon its surface. As a matter of fact, using the d'Alembert's equation, which reads in the present case of spherical symmetry

∂ 2 ∂t 2 -c 2 ∂ 2 ∂r 2
r p(r, t) = 0, (

and the boundary conditions (no singularity of p a (0, t) and continuity p a (R, t) = p w (R, t)), one can in principle derive the air's motion, provided of course that one knows p w (R, t). It is our aim in the present section III to determine p w (R, t) (subsection III A) and to use the result to illustrate the resonance phenomenon (subsection III B).

A. Determining pw(R, t)

Let us denote the common radial displacement of air and water at the bubble's surface by ξ(t):

u a (R, t) = u w (R, t) = ξ(t). (3.2)
In the water, the most general solution of the d'Alembert's equation ((3.1) with c = c w and p = p w ) is the sum of an incoming and an outgoing waves, which can read

p w (r, t) = p in (r, t) + p out (r, t), (3.3a)
where

p in (r, t) = 1 r g t + r -R c w , p out (r, t) = 1 r f t - r -R c w , (3.3b)
with g and f any functions of time. Using Newton's law (2.1b) at r = R, one gets

ρ w ξ = -(grad p w ) r=R = 1 R 2 1 - R c w ∂ ∂t g(t) + 1 + R c w ∂ ∂t f (t) . (3.4)
(i ). To begin with, let us suppose that there is no pressure wave coming from the infinite. Function g should then be taken equal to zero. Solving equation (3.4) for f is straightforward and one gets symbolically [START_REF] Templin | [END_REF]12

f (t) = ρ w R 2 1 + R c w ∂ ∂t (-1) ξ(t) p w (r, t) = p rad (r, t), (3.5a)
where

p rad (r, t) = ρ w R 2 r 1 + R c w ∂ ∂t (-1) ξ t - r -R c w (3.5b)
is the outgoing spherical pressure wave radiated by the pulsating bubble. Not surprisingly, p rad is determined by ξ and its time-derivatives.

(ii ). Let us now consider another situation. We assume that there is a pressure wave p in coming from the infinite (function g is thus nonzero in equations (3.3)-(3.4)), but that the bubble's surface is motionless (function ξ(t) is thus zero ∀t). Solving equation (3.4) for f gives

f (t) = -1 + R c w ∂ ∂t (-1) 1 - R c w ∂ ∂t g(t)
p w (r, t) = p refl (r, t), (3.6a)

where

p refl (r, t) = - R r 1 + R c w ∂ ∂t (-1) 1 - R c w ∂ ∂t p in R, t - r -R c w (3.6b)
is the reflected wave on a rigid spherical obstacle with radius R.

(iii ). Combining the above situations (i ) and (ii ), we finally consider the general case of a pressure wave coming from the infinite onto a pulsating bubble. Solving equation (3.4) for f , we derive the "superposition theorem"

p out (r, t) = p rad (r, t) + p refl (r, t), (3.7)
which is by the way a quite intuitive result.

To conclude this subsection's study, we can now propose a compact expression of the full pressure at r = R. Using (3.5b) and (3.6b) to determine p out (r, t), we get, all calculations carried out,

p w (R, t) = 1 + R c w ∂ ∂t (-1) ρ w R ξ(t) + 2R c w ∂p in (R, t) ∂t . (3.8)
We shall use the above (exact) result in the next section subsection III B and in section IV.

B. Back to the resonances

Let us now consider again the radial eigenmodes of the {air-infinite water} system. Solution (2.4b) can be regarded as the superposition of two monochromatic spherical waves: an incoming wave with complex amplitude p in (r) = -W 2ir e -i(qw(r-R)+ϕ) and an outgoing wave p out (r) = W 2ir e i(qw(r-R)+ϕ) . At r = R, the outgoing wave is thus dephased by an angle π + 2ϕ with respect to the incoming wave:

p out (R) = -e 2iϕ p in (R).
(3.9a) From (3.6b) and the definition (2.5b) of x w , one also derives the dephasing between the reflected wave and the incoming wave:

p refl (R) = - 1 + ix w 1 -ix w p in (R) = -e 2iθ p in (R), (3.9b)
with tan θ = x w . By the way, one can check from expression (2.7) that phase angles ϕ and θ do coincide at the limit χ a → 0 (then x a → 0 and tan ϕ → x w ), as expected for the reflection onto a rigid sphere. Using (3.5b) at r = R and (2.5b) again, one derives the expression of the radiated wave:

p rad (R) = W R sin ϕ -x w cos ϕ 1 -ix w = (e 2iθ -e 2iϕ ) p in (R).
(3.9c)

Adding expressions (3.9b) and (3.9c) gives (3.9a), as expected from the "superposition theorem" (3.7). We have displayed the complex amplitudes p in (R), p out (R), p refl (R) and p rad (R) in figure 3. At a resonance, angle ϕ is equal to π 2 , and then the incoming and outgoing pressure waves are in phase at r = R, so that the bubble's surface corresponds to an antinode of the total pressure wave p w (r). Observe that the incoming wave and the outgoing wave are in phase at r = R if angle ϕ is equal to π 2 . In this case, there is an antinode of the total pressure wave p w (r) at the bubble's surface, which is referred to as a "resonance" in the present paper.

One can add a further remark about the resonances. Let us suppose that, instead of a monochromatic wave, a wave packet is coming from r = ∞. This wave packet can read

p in (r, t) = - 1 2ir dω W(ω) e -i[qw(ω)(r-R)+ϕ(ω)] e -iωt ,
(3.10a) while the corresponding outgoing wave packet is

p out (r, t) = 1 2ir dω W(ω) e i[qw(ω)(r-R)+ϕ(ω)] e -iωt .
(3.10b) Let us suppose, for the sake of simplicity, that the spectral density W(ω) is real. Using the stationary phase condition to determine the position of the wave packet's center, we get

- dq w dω (r -R) - dϕ dω -t = 0 r -R = -c w t + dϕ dω (3.11a)
for the incoming packet, which consequently reaches the bubble's surface at time -dϕ dω , and

dq w dω (r -R) + dϕ dω -t = 0 r -R = c w t - dϕ dω (3.11b)
for the outgoing packet, that leaves the bubble's surface at time dϕ dω . Thus, the wave packet "spends" a time of order ∆t = 2 dϕ dω (> 0) in the bubble. This time is negligible out of resonance where the curve ϕ(ω) (≃ 0 or π) is flat. It is substantial at resonance where ϕ(ω) jumps from 0 to π within an interval ∆ω ≃ χw χa cw R (see (2.11)): thus we have ∆t ≃ 2π ∆ω ≃ χa χw 2R cw , much longer than 2R cw (the corresponding travel time in absence of bubble), which characterizes the resonance as well.

IV. THE MINNAERT BUBBLE

Our aim in this section is to implement the derivation scheduled in introduction and show that, at the isentropic approximation, the Minnaert bubble can be regarded as a one-degree of freedom harmonic oscillator, with a mass, a stiffness, a radiative viscous damping rate and an external forcing term. We begin with the calculations of the total kinetic and (elastic) potential energies of the air bubble. We then derive a motion equation for the "naked" bubble (4.6a), exact but hardly utilizable. This equation is turned, through an unexpectedly exact transformation, in a motion equation for the "dressed" bubble (4.7), readily utilizable.

As suggested by approximated formulas (2.12a,b), one can describe the dynamics of the Minnaert bubble in the following way. Let R+ξ(t) be the bubble's radius at time t. The radial displacement of air simply reads

u a (r, t) = ξ(t) λ a ( r), where λ a ( r) = r R (r R). (4.1)
The above expression of the air displacement is referred to as the "breathing approximation" throughout this paper. As a consequence of (4.1), the total kinetic energy of the air inside the bubble is According to (2.1a) and (4.1), the (homogeneous) extrapressure inside the bubble is

E k = 1 2 bubble ρ a d 3 r ∂ u a ∂t 2 = 1 2 m ξ2 , ( 4
p a = - 1 χ a div u a = - 3ξ χ a R . (4.3)
The classical potential energy of the air inside the bubble thus reads

E p = - ξ 0 4πR 2 dξ ′ p a (ξ ′ ) = 1 2 χ a p 2 a × 4 3 πR 3 = 1 2 Kξ 2 , (4.4a) with K = bubble d 3 r (div λ a ) 2 χ a = 12πR χ a . (4.4b)
We have introduced the mass m and the stiffness K. We should now establish the motion equation ruling the dynamical variable ξ. In order to do it, let us come back to Newton's law for the air; multiplying both hand sides of (2.1b) by λ a , and integrating over the bubble's volume, we get bubble

ρ a d 3 r λ a . ∂ 2 u a ∂t 2 = - bubble d 3 r λ a . -→ grad p = bubble d 3 r (p div λ a -div(p λ a )). (4.5a)
Then, using (2.1a), (4.1) and the Ostrogradsky's theorem, and allowing for definitions (4.2b) and (4.4b), we obtain

m ξ + Kξ = - bubble p d 2 r = -4πR 2 p(r = R, t). (4.5b)
The above motion equation is that of a harmonic oscillator with mass m and stiffness K, submitted to the external (pressure) force 4πR 2 p(r = R, t). Observe that m involves the mere air's inertia (see (4.2b)): equation (4.5b) rules the time-evolution of what could be called the "naked" bubble. Using expression (3.8) of p w (R, t) (= p(R, t)), we get

m ξ + Kξ = -M 1 + R c w ∂ ∂t (-1) ξ(t) + 2 ρ w c w ∂p in (R, t) ∂t , (4.6a)
with 1) ξ term in the right-hand side of equation (4.6a) much looks like an Abraham-Lorentz force, as that undergone by accelerated charged particles in electrodynamics [START_REF] Jackson | Classical Electrodynamics[END_REF] .

M = ρ w 4πR 3 . (4.6b) It is noteworthy that the -M 1 + R cw ∂ ∂t (-
Although exact, equation (4.6a) is nevertheless of poor help to describe the dynamics of the bubble because, even for a free motion, i.e. if no external incoming pressure wave is applied (p in (R, t) = 0), its right-hand side is nonzero, and even large compared to m ξ as we shall see. In the aim of finding a motion equation the right-hand side of which effectively vanishes when p in (R, t) = 0, we can transform equation (4.6a) in the following way. First we apply 1 + R cw ∂ ∂t to both hand sides; then we remove the remaining mass-renormalization term -M ξ and incorporate it with the left-hand side. We get in fine

(m + M ) ξ + RK c w ξ + Kξ + R c w m ... ξ = - 2M ρ w c w ∂p in (R, t) ∂t . (4.7)
The above equation (4.7) is, as far as we know, unpublished. It is that of a harmonic oscillator with mass m + M , stiffness K, and a twofold "friction" mechanism associated with the odd-order time-derivatives of ξ (namely ξ and ... ξ ). This oscillator, that we shall henceforth call the "dressed" bubble, undergoes an external force if -and only if, contrary to the naked bubble -a wave pressure comes from the infinite (see the right-hand side of (4.7)). Observe by the way that this external force term is proportional to the time-derivative of p in (R, t). As can be checked by comparing expressions (4.2b) and (4.6b), the "dressed" mass m+M of the bubble is roughly thousandfold its "naked" mass m: in the acoustic bubble case, the renormalization of the mass due to water's inertia should by no means be regarded as a small effect or treated perturbatively as could be done in the electrodynamics domain. On the other hand, the "friction" terms are in fact radiative terms: they take the emission of spherical acoustic waves by the dressed bubble into account, and they naturally vanish in the incompressible water limit (i.e. for c w → ∞), in which case equation (4.7) reduces to

(m + M ) ξ + Kξ = 0. ( 4.8)
According to the above equation (4.8), the free angular frequency reads

ω 0 = K m + M ≃ K M 1 - m M = ω M 1 - 1 10
ρ a ρ w ≃ ω M , (4.9)
which occurs to coincide with expression (2.9b) of the singular resonance's angular frequency of the {air-water} system (note that this coincidence does not hold for higher orders in the ratio ρ a /ρ w ).

Taking again water compressibility into account, it is easy to determine, at this angular frequency ω 0 , the ratio of the radiative terms in the left-hand side of (4.7). Since ... ξ ≃ -ω 2 0 ξ, and using (4.9), we have

R cw m ... ξ RK cw ξ ≃ mω 2 0 K = m m + M ≪ 1. (4.10)
Dividing equation (4.7) by m + M and neglecting the ... ξ term, we get

ξ + Γ ξ + ω 2 0 ξ = - 2 ρ w c w M M + m ∂p in (R, t) ∂t , (4.11a)
with ω 0 given by (4.9) and

Γ = RK (m + M )c w = Rω 2 0 c w ≪ ω 0 . (4.11b)
Equation (4.11a) describes the motion of a (forced) harmonic oscillator with an effective weak radiative damping rate Γ . Due to the damping, the resonance angular frequency is still (roughly) ω 0 , but its linewidth is finite. In this respect, it is very instructive to calculate the relative linewidth Γ ω0 of this resonance. Owing to (2.9c) and (4.11b), we get

Γ ω 0 = Rω 0 c w ≃ Rω M c w = 3χ w χ a ∼ 6.5 10 -3 (4.12a)
at P 0 = 1 bar. Now, if we consider the relative linewidth of the singular peak in figure 2, we find, using (2.9c) again and (2.11b),

δ x a = 3 ρ a ρ w c a c w c a Rω 0 ≃ 3 χ a ρ w c w Rω M = 3χ w χ a . (4.12b)
Results (4.9) and (4.12a,b) fairly show that the Minnaert model sketched in this section IV is relevant to describe the pulsating bubble physics in the "low frequency" domain, i.e. the frequency range where condition x a ≪ 1 is fulfilled. Furthermore it is readily checked that, within this low frequency domain, expressions (2.7) of sin ϕ and cos ϕ simplify in

sin ϕ ≃ Γ ω (ω 2 0 -ω 2 ) 2 + Γ 2 ω 2 , cos ϕ ≃ ω 2 0 -ω 2 (ω 2 0 -ω 2 ) 2 + Γ 2 ω 2 , (4.13)
which are well known from the forced harmonic oscillator basic theory: in a monochromatic regime with angular frequency ω, angle ϕ is thus the dephasing between the complex amplitudes of ξ and of the forcing term in the right-hand side of (4.11a).

We would add a further remark. In experimental situations, it is unusual that an air bubble is excited by means of a spherical (and moreover concentric) acoustic wave. In practise, a monochromatic plane wave of the type p(z) = p 0 e iqwz is more likely to be used. Now, this plane wave can be expanded on the spherical waves basis 10 : p 0 e iqwz = p 0 e iqwr cos θ = p 0 ∞ ℓ=0 i ℓ (2ℓ + 1)j ℓ (q w r)P ℓ (cos θ), (4.14) where the j ℓ are the spherical Bessel's functions and the P ℓ are the Legendre's polynomials. In fact, in the "low frequency" domain (as defined above), we have q w R ≪ 1. Now, only j 0 (ρ) = sin ρ ρ is non negligible for ρ ≪ 1. Consequently, only the ℓ = 0 wave component in expansion (4.14) will significantly affect the bubble. Said in other words: as far as the wavelength is large compared with its radius R, the bubble makes no difference between a plane wave and a spherical wave. As a conclusion, the motion equation of an air bubble located at the origin of coordinates and undergoing the monochromatic pressure wave p 0 e iqwz will read (see (4.11a))

(ω 2 0 -ω 2 -iωΓ ) ξ ≃ - p 0 ρ w R e -iqwR . (4.15a)
Expanding the above e -iqwR term up to first order in q w R, we finally find, using (4.11b),

ξ ≃ - p 0 ρ w R 1 ω 2 0 -ω 2 -i Rω 3 cw , (4.15b)
with the ω 3 radiative damping term usually found in literature [START_REF] Leighton | The Acoustic Bubble[END_REF]5 . Observe that the discrepancy between the radiative damping rates Γ = Rω 2 0 cw and Γ (ω) = Rω 2 cw is in fact but apparent: the viscous constant rate Γ is found whenever the incident pressure wave is appreciated at r = R (see (4.11a) or (4.15a)); the ω-dependent rate Γ (ω) is found whenever this incident pressure is appreciated at r = 0 (see (4.15b) and, more generally, the scattering functions used in the literature).

How is the above equation modified when dissipation is taken into account ? As already mentioned in subsection II B about the general motion, an exact answer to this question would imply an exact solution of equations (2.18) (under the circumstances, in the x a ≪ 1 limit and in monochromatic regime). This delicate enterprise, important though it may be in other respects, is outside the scope of the present study. To manage it, Prosperetti proposes a detailed discussion, 5 using the following simplifications.

(i ). The water is assumed to be incompressible as concerns viscosity effects : the Navier-Stokes equation (2.18b) thus reduces to (2.1b) in the bulk; at the interface (r = R), the viscous stress results in an additional extrapressure p vis = 4η w ξ R which in fine finds expression in a viscous damping rate Γ vis = 4ηw ρwR 2 . (ii ). Due to the water's high heat capacity, the extratemperature θ(r, t) is neglected for r R: the water is regarded as a heat reservoir with temperature T 0 .

(iii ). Any viscosity of the gas is neglected : equation (2.1b) is available for r < R.

With the above approximations, it can be shown that equation (2.18a) can be heuristicly simplified in p(r) = -1 χ a (ω) div u(r), with χ a (ω) = χ a γ κ(ω) (κ(ω) is the socalled polytropic exponent). Accordingly, the resonance frequency ω 0 is now frequency-dependent and reads ω 2 0 (ω) = κ(ω) γ ω 2 M : the Minnaert bubble is, strictly speaking, no longer a harmonic oscillator in the sense that its stiffness K(ω) = 12πR χ a (ω) is frequency-dependent (observe that its mass M = 4πρ w R 3 is unaffected by heat transfers in the gas). Moreover, thermal effects yield an additional frequency-dependent thermal damping, so that the left-hand side of equation (4.15a) should read in fine ω 2 0 (ω) -ω 2 -iωΓ eff (ω), with Γ eff (ω) = Γ + Γ vis + Γ th (ω). For a bubble with radius R = 1 mm, one has Γ = 260 s -1 , Γ vis = 7.2 10 -2 s -1 , Γ th (ω M ) = 540 s -1 : viscosity is fairly negligible, and the damping is dominated by the thermal process. Nevertheless, in the next section, we shall ignore dissipation again for the sake of simplicity.

V. THE ACOUSTIC RADIATIVE PICTURE APPROACH

As discussed in the previous section, 99.9 % of the dressed bubble's mass is due to water's inertia. The dressed bubble is consequently, at least partially, "made of water". It is thus tempting to raise the naive question: what part of the water is devoted to the dressing of the bubble, and what part to the propagation of sound? Our aim in this last section is to provide an answer to this puzzling question. We propose an alternative description of the Minnaert bubble, including a possible departure from the breathing approximation for the air motion, as well as a physical interpretation of the mass-renormalization term. To begin with, we introduce an ad hoc splitting of the displacement field u( r, t).

(i ). If the breathing approximation was exact, the air displacement would read like (4.1). In order to take a possible departure from this behaviour into account, let us introduce a further degree of freedom for the air and set

u a ( r, t) = ξ(t) λ a ( r) + u ar ( r, t) (r R) (5.1)
( u ar thus vanishes at the breathing approximation).

(ii ). If the water was incompressible, its displacement would read u wi ( r, t) = ξ(t)λ w ( r), where λ w ( r) = R 2 r 2 e r (r R).

(5.2a)

In order to take water compressibility into account, let us write the water displacement in the form u w ( r, t) = u wi ( r, t) + u wr ( r, t) (r R) (5.2b) ( u wr thus vanishes in the incompressible water limit).

The displacement field u r ( u ar for r R, u wr for r R) will be referred to as the "radiative" displacement field. It is noteworthy that the above splitting of the water displacement in an (instantaneous) "elastostatic" field ξ(t) λ w ( r) and a "radiative" field u wr ( r, t) is the acoustic transposition of the well known splitting of the total electric field in an electrostatic part E e = -∇V and a radiative part E r = -∂ A ∂t in the electromagnetism domain. Let us recall that the electrostatic field is instantaneous in the Coulomb gauge (div A = 0), also called "radiation gauge". Hence our coining the term "radiative". Of course, our splitting of the air displacement meets no electromagnetic equivalent, since the interior of a charged particle does not make sense. Now, let us consider the lagrangian describing the radial motion of the {air-water} system:

L = ∞ 0 4πr 2 dr 1 2 ρ ∂ u ∂t 2 - 1 2χ (div u) 2 , (5.3)
with (ρ, χ) = (ρ a , χ a ) for r < R, (ρ w , χ w ) for r > R.

With the splittings (5.1) and (5.2b), the above lagrangian (5.3) also reads exactly

L = L b + L r + L br , (5.4)
with

L b = 1 2 (m + M ) ξ2 - 1 2 Kξ 2 , (5.5a) L r = ∞ 0 4πr 2 dr 1 2 ρ ∂ u r ∂t 2 - 1 2χ (div u r ) 2 ,(5.5b) L br = ξ ∞ 0 4πr 2 drρ λ. ∂ u r ∂t . (5.5c)
Let us comment these three terms.

-Lagrangian L b describes the free motion of the Minnaert bubble as a harmonic oscillator with mass m+M and stiffness K. We should nevertheless emphasize the following point: if M is fairly well the same mass ρ w 4πR 3 introduced in equations (4.6), it no longer corresponds here to an Abraham-Lorentz coefficient as in section IV, but to the value of the integral

M = ∞ R 4πr 2 dr ρ w λ 2 w (r).
(5.6a)

In other words, 1 2 M ξ2 is the total kinetic energy of the water surrounding the bubble in the incompressible water limit. Hence the very intuitive physical meaning of the mass renormalization put forward by Minnaert himself 1 . In this respect, it should be remarked that, in the case of a finite extent of water considered in subsection II A, integral (5.6a) should be substituted by

M (L) = L R 4πr 2 dr ρ w λ 2 w (r) = 1 - R L M (< M ), (5.6b)
which gives for the Minnaert harmonic oscillator the L-dependent angular frequency

ω M (L) = K M (L) = L L -R ω M (> ω M ) (5.7)
to be compared with expression (2.24).

-Lagrangian L r describes the free radial motions of air (interval [0, R]) or water (interval [R, ∞[), these motions being uncoupled de facto from each other by the rigid boundary condition at r = R: u ar (R, t) = u wr (R, t) = 0, ∀t.

-Lagrangian L br couples the Minnaert bubble and the radiative field. Observe that, due to the timederivatives ξ and ∂ ur ∂t , this coupling is of the inertial type.

From lagrangian L, one can derive the equations of motion, as well for the Minnaert bubble as for the radiative sound field.

-For the Minnaert bubble, we have

(m + M ) ξ + Kξ = - ∞ 0 4πr 2 dr ρ λ. ∂ 2 u r ∂t 2 .
(5.8)

The fact that the right-hand side of the above equation is an integral is easily understandable: according to splittings (5.1) and (5.2b), the fluid motion is "shared out" among the Minnaert motion and the radiative motion; consequently, whereas the air bubble is strictly confined to the 0 < r < R domain, the Minnaert bubble stretches to infinity.

-For the radiative sound field, we have

∂ 2 u r ∂t 2 -c 2 -→ grad(div u r ) = -ξ λ.
(5.9)

As expected from the electromagnetic analogy with the radiation of a moving charge, the source term of the excitation of the radiative acoustic field u r is proportional to the second time-derivative of dynamical variable ξ. Observe too that, in perfect symmetry with the above-discussed excitation of the bubble by the sound field, this source-term (-ξ λ) is not localized at the air-water boundary but spreads inside the air and dies out in the water like λ w ( r), i.e. like 1/r 2 .

We would finish this discussion with a further remark. Let us come back to expression (3.5b) of the water extrapressure p rad (r, t). If water was incompressible, implying c w = ∞, equation (3.5b) would simply read p w (r, t) = p wi (r, t) = ρ w R 2 r ξ(t).

(5.10a)

It is then tempting to parallel the splitting (5.2b) of the water displacement, and introduce the radiative extrapressure p wr defined by p w (r, t) = p wi (r, t) + p wr (r, t).

(5.10b)

We have then, allowing for (5.2a),

ρ w ∂ 2 u wi ∂t 2 = ρ w ξ R 2
r 2 e r = --→ grad p wi , (5.10c) and consequently (Newton's law):

ρ w ∂ 2 u wr ∂t 2 = --→ grad p wr . (5.10d)
In this respect, it is remarkable that Newton's law (2.1b) applies separately for each couple ( u wi , p wi ), ( u wr , p wr ). This is due to the fact that extrapressure p wi is purely inertial, i.e. archimedean, the water acceleration ξ λ w standing for a (local) pseudo gravity.

VI. CONCLUSION

The {air bubble-water} system has a continuum of radial eigenmodes: for any given angular frequency ω, one can build (subsection II A) a spherical standing wave satisfying the equation of motion and the boundary conditions. For a discrete series of values of ω, this standing wave has, in the water side, a pressure antinode at r = R (bubble's surface). Throughout this paper, these particular modes are referred to as resonances. They can all the same be understood as the resonances of an air-water spherical Pérot-Fabry cavity. Due to this spherical symmetry, the fundamental resonance's angular frequency ω 0 is spectacularly lowered with respect to what would be expected in a 1-D plane geometry, and is found to coincide with the so-called Minnaert angular frequency. We have also considered (subsection II C) the case of an airbubble with a finite extent of surrounding water, and discussed the conditions for the infinite water approximation to be available. Within the framework of this approximation, we have derived (subsection III A) an exact expression for the pressure exerted by the water onto the bubble's surface, and interpreted the resonances from the point of view of the water (subsection III B). At low enough frequencies, the acoustic wavelength (in the air) is large compared to the bubble's radius, so that the motion of the air can be approximated by a homogeneous inflation/deflation of the bubble. At this breathing approximation, the bubble's dynamics is described by a single dynamical variable ξ, which is shown to be ruled by a harmonic oscillator motion equation. We have shown how the "naked" bubble can be "dressed" by the pressure forces exerted on its surface by the outer water, and we finally conclude that the Minnaert bubble is a har-monic oscillator slightly different from the one generally considered in the literature (section IV). In order to interpret the mass renormalization term that appears with the dressing of the naked bubble, we have outlined (section V) an alternative picture, which consists in splitting the fluid (air or water) displacements in two components, as currently done in electrodynamics when splitting the full electric field in an electrostatic (coulombean) and a radiative parts. Besides, this "radiative" picture can be readily extended to the study of many-bubble clouds.

FIG. 2 :

 2 FIG.2: sin ϕ as a function of xa = qaR. In addition to regular peaks for xa ≃ 3π 2 , 5π 2 , . . . , (2n+1)π 2

  the solution ω p of equation (2.20) lies within the interval [(p -1) cwπ L-R , p cwπ L-R ], and thus, for L → ∞, the continuum is recovered. Moreover, the pressure wave in water associated with eigenmode p reads (see (2.4b)) p w (r) = W r sin[ ωp cw (r-R)+ϕ(ω p )]. It consequently displays an antinode (pressure loop) at r = R provided that ϕ(ω p ) = π 2 , in which case cot( ωp cw (L -R)) = 0 and equation (2.21) reduces to (2.8). If the maximum distance between two consecutive eigen angular frequencies is small compared to the linewidth of function sin ϕ displayed in equations (2.11), namely if one has

  remains much smaller than unity, the expansion in equation (2.9a) still holds and we calculate a L-dependent Minnaert angular frequency

FIG. 3 :

 3 FIG.3: Complex amplitudes of the incoming (p in ) and outgoing (p out , p refl and p rad ) waves at r = R, illustrating the "superposition theorem" (3.7) in the monochromatic regime. Observe that the incoming wave and the outgoing wave are in phase at r = R if angle ϕ is equal to π 2 . In this case, there is an antinode of the total pressure wave p w (r) at the bubble's surface, which is referred to as a "resonance" in the present paper.

.2a) with m = bubble ρ a d 3 r λ 2 a = ρ a 4 5 πR 3 .

 23 (4.2b)

11

 11

  1) can be regarded as the symbolic summation of the infinite series 1 -R cw ∂ ∂t + `R cw ´2 ∂ 2 ∂t 2 + . . . .

  .22) then we are sure to find all the resonances foreshadowed by equation (2.8) among the set of eigen angular frequencies satisfying (2.20). It is noteworthy that the above condition (2.22) is extremely severe. With χ w ≃ 10 -10 Pa -1


Bruno Denet
email: bruno.denet@irphe.univ-mrs.fr

Sivashinsky equation in a rectangular domain

The (Michelson) Sivashinsky equation of premixed flames is studied in a rectangular domain in two dimensions. A huge number of 2D stationary solutions are trivially obtained by addition of two 1D solutions. With Neumann boundary conditions, it is shown numerically that adding two stable 1D solutions leads to a 2D stable solution. This type of solution is shown to play an important role in the dynamics of the equation with additive noise.

I. INTRODUCTION

The Sivashinsky equation [1] (or Michelson Sivashinsky equation depending on the authors) is a non linear equation which describes the time evolution of premixed flames. Because of the jump of temperature (and thus of density) across the flame, a plane flame front is submitted to a hydrodynamic instability called the Darrieus-Landau instability. The conservation of normal mass flux and tangential velocity across the front leads to a deflection of streamlines which is the main cause of this instability. A more detailed description of this instability can be found in the book of Williams [START_REF] Williams | Combustion Theory[END_REF] (see also, in the approximation of potential flow in the burnt gases, the elementary electrostatic explanation in [START_REF] Denet | [END_REF], where essentially the flame is described as a surface with a uniform charge density). At small scales, the instability is damped by diffusive effects: the local front propagation velocity is modified by a term proportional to curvature, the coefficient ahead of the curvature term is called the Markstein length. A geometrical non linear term, which limits ultimately the growth of the instability, is caused by the normal propagation of the flame. The Sivashinsky equation, obtained as a development in powers of a gas expansion parameter, i.e. for a small jump of temperature, or equivalently for a flow almost potential in the burnt gases, represents a balance between the evolution due to these three effects, Darrieus Landau instability, stabilization by curvature and normal propagation of the flame.

The qualitative agreement between the Sivashinsky equation and direct numerical simulations, generally performed with periodic boundary conditions has been excellent even with large gas expansion [4], and also when gravity is included [5]. It has been shown in a classic paper of the field [6] (following [7], where the pole decomposition was introduced) that the 1D solution of the Sivashinsky equation in the absence of noise was attracted for large times toward stationary solutions, with poles aligned in the complex plane, called coalescent solutions. It was shown analytically in [8] [9] that each solution, with a given number of poles, is linearly stable in a given interval for the control parameter (either the domain width or more often the Markstein length with a domain width fixed to 2π).

In a recent paper [10] (hereafter called I) the present author has been interested in the behavior of the Sivashinsky equation in 1D, but with Neumann boundary conditions (zero slope at each end of the domain), a situation which, although more realistic than periodic boundary conditions, had not attracted much interest over the years. Actually, periodic boundary conditions lead to a symmetry which is not present in the case of a flame in a tube, i.e. every lateral translation of a given solution is also a solution. Presented in a different way, a perturbation on the flame can grow, reach the cusp (the very curved part of the front, pointing toward the burnt gases) and then decay, but after having caused a global translation of the original solution. This is not possible with Neumann boundary conditions, but it was supposed that this difference with periodic boundary conditions was unimportant. The surprise was however that stable stationary solutions in the Neumann case involved a number of solutions with two cusps (and the corresponding poles) at each end of the domain, called bicoalescent solutions. This type of solution of the Sivashinsky equation was already introduced in [11], although this last article did not obtain those which are stable with Neumann boundary conditions. The author would like to mention here two articles which he did not cite in I , namely [12], where some bicoalescent solutions with Neumann boundary conditions were first obtained, and [13] where bicoalescent solutions were obtained in direct numerical simulations. In this last paper, one solution was not completely stationary, because of the effect of noise, but another solution was actually almost stationary. Of course the computer time needed for such a simulation is probably one hundred times more than the equivalent Sivashinsky equation simulation, with all sorts of possible sources of noise, so obtaining really stationary bicoalescent solutions in this case is a challenging task.

Coming back to I, we can summarize the 1D results of this paper in the following way:

1. Bicoalescent solutions were obtained, stable with Neumann boundary conditions. Simulations performed without noise tend to these solutions.

2. The new solutions led to a bifurcation diagram with a large number of stationary solutions, where particularly the number of solutions multiply when the Markstein length, presented above, which controls the stabilizing influence of the curvature term, decreases 3. The bicoalescent solutions play a major role in the dynamics of the equation with additive noise. In the case of moderate white noise, the dynamics is controlled by jumps between different bicoalescent solutions.

In the present paper, we shall be interested in the Sivashinsky equation with Neumann boundary conditions, but in two dimensions in a rectangular domain. Another nice property of the equation (apart from the pole decomposition in 1D) is that 2D solutions can be formed by the simple addition of two 1D solutions, one for each coordinate [START_REF] Joulin | Image des Mathématiques Modélisation de la combustion[END_REF]. The exact counterpart of I will be obtained:

1. Sums of two bicoalescent solutions are stable in 2D with Neumann boundary conditions. The time evolution of the equation without noise tends toward these solutions 2. With sums of a large number of 1D solutions, a really huge number of 2D solutions can be obtained.

3. The sums of bicoalescent solutions play also a major role in the dynamics in two dimensions in the presence of noise.

II. SOLUTIONS IN ELONGATED DOMAINS

The Sivashinsky equation in one dimension can be written as

φ t + 1 2 φ 2 x = νφ xx + I (φ) (1)
where φ (x) is the vertical position of the front. The Landau operator I (φ) corresponds to a multiplication by |k| in Fourier space, where k is the wavevector, and physically to the destabilizing influence of gas expansion on the flame front (known as the Darrieus-Landau instability, and described in the introduction). ν is the only parameter of the equation (the Markstein length) and controls the stabilizing influence of curvature. The linear dispersion relation giving the growth rate σ versus the wavevector is, including the two effects:

σ = |k| -νk 2 (2)
As usual with Sivashinsky-type equations, the only non linear term added to the equation is 1 2 φ 2 x . In the flame front case, this term is purely geometrical : the flame propagates in the direction of its normal, a projection on the vertical (y) direction gives the factor cos (θ) = 1/ 1 + φ 2

x , where θ is the angle between the normal and the vertical direction, then a development valid for small slopes of the front leads to the term 1 2 φ 2

x . The Sivashinsky equation is typically solved numerically on [0, 2π] with periodic boundary conditions. In I it has also been solved on [0, 2π] with only symmetric modes, which corresponds to homogeneous Neumann boundary conditions on [0, π] (zero slope on both ends of the domain).

The two dimensional version of the Sivashinsky equation is

φ t + 1 2 (∇φ) 2 = ν∆φ + I (φ) (3)
where the Landau operator I (φ) corresponds now to a multiplication by k 2 x + k 2 y in Fourier space, k x and k y being the wavevectors in the x and y directions. All dynamical calculations, are performed by Fourier pseudo-spectral methods (i.e. the non linear term is calculated in physical space and not by a convolution product in Fourier space). The method used is first order in time and semi-implicit (implicit on the linear terms of the equation,

φ = 2ν N n=1 ln sin x -z n (t) 2 + ln sin x -z * n (t) 2 ( 4
)
where N is the number of poles z n (t) in the complex plane. Actually the poles appear in complex conjugate pairs, and the asterisk in Equation 4 denotes the complex conjugate. In all the paper, the number of poles will also mean number of poles with a positive imaginary part. The pole decomposition transforms the solution of the Sivashinsky equation into the solution of a dynamical system for the locations of the poles. In the case of stationary solutions, the locations of the poles are obtained by solving a non linear system:

-ν 2N l=1,l =n cot z n -z l 2 -isgn [Im (z n )] = 0 n = 1, • • • , N (5)
where Im (z n ) denotes the imaginary part and sgn is the signum function. This non linear system is solved by a variant of Newton method.

With periodic boundary condition, the usual result is that in the window 2n -1≤ 1/ν≤ 2n+ 1, n = 1, 2, • • • there exists n different monocoalescent stationary solutions (all the poles have the same real part), with 1 to n poles, and the solution with the maximum number of poles n is asymptotically stable. For a particular value of 1/ν, the number n(ν) such that 2n -1≤ 1/ν≤ 2n + 1 is called the optimal number of poles.

With Neumann boundary conditions, in each of the intervals [2n -1, 2n + 1] of the parameter 1/ν, not only one asymptotically stable solution, but n + 1, of the form (l, n -l)

with l = 0, 1,• • • , n where l poles coalesce at x = 0 and l -n coalesce at x = π, were obtained in I. (The bicoalescent type of solutions have been recently introduced in [11]). In Figure 1 is shown a bifurcation diagram with all the possible stable stationary solutions (plotted only, contrary to I, in the domain where they are stable) versus 1/ν. What is actually plotted is the amplitude ∆φ (maximum minus minimum of φ) versus 1/ν. As can be seen, when the optimal number of poles increases with 1/ν, the number of stable stationary bicoalescent solutions is also increasing. The stability of these solutions is not proved analytically, nor by a numerical study of the linearized problem, we use only numerical simulations of the Sivashinsky equation, with the different bicoalescent solutions plus some small perturbations as initial conditions, and the solution returns toward the unperturbed solution.

In a square domain [0, 2π] * [0, 2π], it has been remarked in [START_REF] Joulin | Image des Mathématiques Modélisation de la combustion[END_REF] that if φ 1 (x) and φ 2 (x) are solutions of the 1D Sivashinsky equation ( 1), then φ 1 (x) + φ 2 (y) (we use here

φ 1 ⊕ φ 2
as a notation for this sum, whose amplitude is the sum of the amplitudes of φ 1 and φ 2 ) is a solution of (3) in two dimensions, and that the stationary solution obtained numerically in this case for periodic boundary conditions [START_REF] Michelson | [END_REF] is simply a sum of two monocoalescent 1D solutions. Let us note that, if it is absolutely obvious that sums are solutions of the 2D equation, the stability of these solutions has never been proved analytically, and can only be inferred from a small number of numerical simulations.

In the case of rectangular domains [0, 2π] * [0, 2b], sums are also solutions of the equation, with φ 2 now solution of the 1D Sivashinsky equation with parameter 1/ν in domain [0, 2b],

which can be obtained by an appropriate rescaling from the solution in [0, 2π] with parameter

1/ν 1 = (1/ν)(b/π).
A particularly simple case is the limit where b is very small, where the only solution with to noise (although less than the pure 1D solutions) it could be possible to observe in direct numerical simulations and experimentally the solutions with the lower amplitude, which are the least sensitive to noise. In experiments, the solutions should also survive heat losses (important in narrow channels) and not be too much perturbed by gravity (i.e. have a large enough Froude number) in order to be observed .

III. SOLUTIONS IN SQUARE DOMAINS

We now turn to stationary solutions in square domains [0, 2π] * [0, 2π] with Neumann boundary conditions. Sums of bicoalescent solutions produce also in this case stable stationary solutions. The purpose of this section is to give details on the consequences of this simple addition property. We show first the different types of solutions obtained by addition of stable bicoalescent solutions in 1D. We insist on the fact that these solutions are linearly stable and give a specific example of the time evolution of one such solution with some small perturbations. Finally two bifurcation diagrams are provided, one is the 2D equivalent of Figure 1 with the stable solutions plotted only in their stable domain. The second contains all the solutions obtained by addition of all the branches found in 1D in I, and as the reader will see, a really huge number of branches are created in this way.

In Figures 3 and4 In Figure 5, we have an example showing the stability of the (3, 2) ⊕ (4, 1) solution. We start from this solution and add an additive white noise to Equation (3) when the time is below 0.5. This white noise is gaussian, of deviation one, and we multiply it by an amplitude a = 0.001. It can be seen that after the noise is stopped, the solution tends exponentially back toward the (3, 2) ⊕ (4, 1) solution. Similar figures would be obtained with the other solutions of Figures 3 and4, except that higher amplitude solutions would need an even lower noise in order not to jump immediately toward a lower amplitude solution.

In Figure 6 is shown the strict 2D equivalent of Figure 1: the bifurcation diagram showing the amplitude versus 1/ν for all the solutions which are linearly stable, only plotted in their domain of stability. For 1/ν < 3 there is only one possibility (1, 0) ⊕ (1, 0). For 3 < 1/ν < 5

we have three branches (from higher to lower amplitudes) (2, 0) ⊕ (2, 0) (1, 1) ⊕ (2, 0) and

(1, 1) ⊕ (1, 1). For the value 1/ν = 10 we have the six solutions of Figures 3 and4, that is from higher to lower amplitudes the (5, 0) ⊕ (5, 0), ( 4 the pole decomposition is not available, but nevertheless a lot of 1D stationary solutions have been obtained [16]. The Kuramoto-Sivashinsky equation shares with the Sivashinsky equation the possibility to create 2D solutions by adding two 1D solutions, so actually in this case we have also a very large number of branches. These rectangular solutions are not as physically relevant in the Kuramoto-Sivashinsky equation case. Contrary to the Sivashinsky equation, where stable stationary solutions are basically as large as possible and are thus rectangular in a rectangular domain, it seems likely that in the Kuramoto-Sivashinsky case, the most interesting solutions would have an hexagonal symmetry (hexagonal cells are also observed for the Sivashinsky equation with stabilizing gravity [17]). Stationary solutions of the Sivashinsky equation with hexagonal symmetry should exist too, and the author conjectures that the order of magnitude of the number of solutions with hexagonal symmetry should be approximately the same as those with rectangular symmetry. Apparently there is no trivial way to construct hexagonal solutions, so unfortunately, until some progress is made, obtaining the hexagonal equivalent of Figure 7 is almost impossible. We have here an example emphasizing the fact that as the smoothing effect (viscosity, curvature, surface tension ...) decreases, we are not able to generate correctly all the simple solutions of a given set of partial differential equations (Sivashinsky and Kuramoto-Sivashinsky equations, Navier Stokes ...) even with the aid of computers.

IV. EVOLUTION WITH NOISE

In the previous section we have shown numerically that the sums of linearly stable 1D bicoalescent solutions lead to linearly stable 2D solutions. However, even a linearly stable solution could have a very small basin of attraction. So in this section, we study the effect of noise on the solutions of the Sivashinsky equation in a square domain, with Neumann boundary conditions. The important solutions will be the solutions that are reasonably resistant to the applied noise.

This noise used here is simply an additive noise, added to the right-hand side of Equation (3). We choose the simplest possible noise, a white noise (in space and time), which is gaussian, has deviation one and is multiplied by an amplitude a. But contrary to Figure 5, this noise will be applied at each time step. We use in all the simulations presented the same parameter 1/ν = 10, the stationary solutions corresponding to this parameter have been presented in the previous section. We recall that in I, for the one dimensional version of the Sivashinsky equation with moderate noise, the evolution was analysed in terms of jumps between the available bicoalescent stationary solutions. We would like to show here that in 2D, the sums of bicoalescent solutions also play an important role in the dynamics.

In Figure 8, starting from an initial condition which is the (4, 1)⊕(4, 1) stationary solution, is plotted the amplitude of the solution versus time, for a noise amplitude a = 0.01, with also straight lines corresponding to the amplitudes of the lowest amplitude stationary solutions, i.e. those of Figure 4. The stationary solutions with higher amplitudes (those of Figure 3) apparently are too sensitive to noise to play any role in the dynamics. It is seen in Figure 8 that because of the noise, the solution departs quickly from the (4, 1) ⊕ (4, 1) solution, and that it seems that, during the time evolution, the solution is close (apart from some violent peaks in the amplitude) to the (3, 2) ⊕ (4, 1) solution for some time, then finally the amplitude decreases again to be near that of the (3, 2) ⊕ (3, 2) solution.

In order to prove that the solution is indeed close to the previously mentioned solutions, because after all very different solutions could have similar amplitudes, we plot in Figure 9 for the same simulation, what we have called the distance between the solution at a given time and the sums of bicoalescent solutions, which is simply the L1 norm of the difference between both solutions. The spatial mean value of all solutions is adjusted here to have the same value. Normally, it is necessary to measure the distance between the solution and all symmetries of a given sum of bicoalescent solutions (i.e. you can interchange the poles at 0 and π in the x and y directions) but for the low amplitude a = 0.01, it has not been necessary, and we plot only the distance from the relevant solutions.

As we start from the (4, 1) ⊕ (4, 1) solution, the distance to this solution is zero initially, and we can see that, although the amplitude seems to indicate that at some time, one is again close to this solution, this is not the case. On the contrary, the solution returns regularly close to the (3, 2) ⊕ (4, 1) solution for times lower than 110, then there is a transition toward something close to the (3, 2) ⊕ (3, 2) solution, the solution departs slightly from this last solution for some time, possibly toward a linearly unstable stationary solution, and returns toward it at the end of the simulation. As Figure 6 remotely looks like energy levels in atomic physics, one could be tempted to interpret the evolution of the two previous figures with a small noise (apparently in 2D the solution is less sensitive to a given amplitude of the white noise compared to 1D simulations) as a sort of deexcitation from the high amplitude level (4, 1)⊕(4, 1) toward first (3, 2)⊕(4, 1), then toward the fundamental level (3, 2)⊕ [START_REF] Denet | [END_REF][START_REF] Williams | Combustion Theory[END_REF].

Indeed, between the sums of bicoalescent solutions, if all are linearly stable, the solutions with the lower amplitude seems to be more resistant to the action of noise.

To better understand the effect of noise, we present now a simulation with a larger noise amplitude a = 0.1 , ten times larger than the previous case (we recall that this noise amplitude should be compared to the laminar flame velocity, which is normalized to 1 in this paper). In Figure 10 is plotted the amplitude versus time, with as before straight lines with the amplitude of the important sums of bicoalescent solutions. The initial condition is also the (4, 1) ⊕ (4, 1) solution. Apparently this last solution is too sensitive to noise to play a meaningful role in the dynamics, although it happens that some peaks in the amplitude could involve solutions not too far from this initial solution. As the distance to this solution is never really small, even in the peaks, we shall not comment further on this solution. On the other hand, it seems that a lot of time is spent with an amplitude close to that of the (3, 2) ⊕ (3, 2) solution (which we have called previously the fundamental level), and perhaps some time with an amplitude close to the (3, 2) ⊕ (4, 1) solution (the first excited level).

In order to see what is really occurring, we now turn to figures of the distance (defined above) to these two solutions versus time (for the same simulation of Figure 10). However, for a higher amplitude, we have to include the four different symmetries of these solutions in the analysis (i.e. for instance (

. In Figure 11 is shown the distance to the four symmetries of the (3, 2) ⊕ (3, 2) solution.

The distance to one of the four symmetries is indeed often small (but not very small for this value of the noise) during the time evolution. Then because of the noise, perturbations are created that lead the amplitude to increase as the perturbation is convected toward one of the cusps, and the solution often comes back toward another symmetry of the fundamental level.

In Figure 12 is shown the distance to the four symmetries of the (3, 2) ⊕ (4, 1) solution (the first excited level) (always for the same simulation). It is seen that the solution is only reasonably close to this type of solution at times close to 50. At other times, minima of the distance are not very small and the solution is often closer to the (3, 2) ⊕ (3, 2) solution. In I, we have presented the evolution of the 1D Sivashinsky equation with a moderate additive noise as a series of jumps between bicoalescent solutions. In 2D the situation is relatively similar, with the sums of bicoalescent solutions playing the same role. However, the noise amplitude necessary to cause jumps seems much higher in 2D, and practically speaking during the previous simulation, only the fundamental (the solution with the lowest amplitude) and first excited levels were obtained. It should also be noted that the degenerescence (the four possible symmetries) of the fundamental level is probably important in the evolution (for instance for 1/ν = 12 the fundamental level would be (3, 3) ⊕ (3, 3), which does not lead to other solutions by symmetry, so that it should be less likely to obtain the fundamental level in this case).

Before closing this section, let us insist on the fact that, if the solution regularly returns toward sums of bicoalescent solutions, the fronts we obtain are not sums for each time. Figure 13, where is plotted a front of the previous simulation for time 120.1 (just before a peak of the amplitude in Figure 10) , should be a clear example of this property. In this Figure ,   the whole domain [0, 2π] * [0, 2π] is plotted as before, but this time as a grayscale figure, white corresponding to the minimum of φ, black to the maximum. Essentially, an oblique perturbation has grown on a front that was previously a sum. This oblique perturbation moves toward each corner of Figure 13, and the amplitude peak corresponds to the moment where the perturbation reaches the corner. Then the solution is attracted again toward a sum of bicoalescent solutions.

To summarize this section on the effect of noise, the fact that all the sums of bicoalescent solutions with the optimal number of poles are linearly stable does not prove that they can be practically observed. On the contrary, the solutions with the larger amplitude have a basin of attraction so small that they can almost never be seen. We have introduced an analogy with atomic physics by calling the bicoalescent solution with the lowest amplitude the fundamental level, other solutions the excited levels. In the examples shown, only the fundamental and first excited levels (and their symmetries) were obtained during the time evolution of the Sivashinsky equation excited by an additive noise. We recall that in I, it was shown in 1D that the evolution with noise was completely different with periodic boundary conditions, where only the largest amplitude monocoalescent solution was linearly stable (even if extremely sensitive to noise). In this case, the solution regularly returns close to the highest amplitude solution. With Neumann boundary conditions, this is just the opposite, the solution prefers to be close to the lowest amplitude, almost symmetric, sum of bicoalescent solutions.

V. CONCLUSION

In this paper, we have used the possibility to create two dimensional rectangular stationary solutions from the addition of two 1D stationary solutions in order to generate a huge number of stationary solutions of the Sivashinsky equation. With Neumann boundary conditions, the addition of two stable 1D bicoalescent solutions leads to stable 2D solutions, which also play a role in the dynamics when an additive noise is added to the equation. However, with noise, only the sums of bicoalescent solutions with the lowest amplitude (which are less sensitive to noise) have a reasonable chance to be observed. More precisely, jumps between different symmetries of the lowest amplitude sum, or between the two sums with the lower amplitude, are obtained in the simulations. Although we have used a white noise in this paper, experiments, submitted to a residual turbulence, should behave in a similar way. In order to have a large enough Froude number for gravity effects to be negligible, flames with a sufficiently large laminar flame velocity would have to be chosen. Presence of an oblique perturbation which has grown on a sum of bicoalescent solutions. This perturbation will reach one corner in the figure, be damped, and the solution will again be close to a sum of bicoalescent solutions.

explicit on 1 2 φ 2 x

 2 ). No particular treatment of aliasing errors is used. The 2D Sivashinsky equation is solved in [0, 2π] * [0, 2b] with only symmetric modes, which corresponds to homogeneous Neumann boundary conditions in the rectangular domain [0, π] * [0, b] . Pole solutions ([6]) of the 1D Sivashinsky equation are solutions of the form:

parameter 1 /

 1 ν in domain [0, 2b] is simply the flat (0) solution φ 2 = 0. As a sum of the previously described bicoalescent solutions in [0, 2π] added to the flat solution in the other direction, we have simply a way to observe the bicoalescent solutions in two dimensions. We have observed numerically (not shown here, the behavior is very similar to the 1D case) for Neumann boundary conditions, that these sums (l, n -l) ⊕ (0) are stable . As an example, for 1/ν = 10 and b = π/10 we show in Figure 2 a perspective view of the three different stationary bicoalescent solutions (5, 0) ⊕ (0), (4, 1) ⊕ (0), (3, 2) ⊕ (0) (from top to bottom). In all the figures, the whole domain [0, 2π] * [0, 2b] is plotted, the solution with Neumann boundary conditions corresponds only to one fourth of the domain [0, π] * [0, b] . We have found it clearer to show the whole domain (contrary to I), because some solutions are very difficult to distinguish if plotted in [0, π] * [0, b]. Although these solutions are very sensitive

  are shown the six stable solutions obtained from the three 1D solutions of Figure2for 1/ν = 10. In Figure3can be seen (in perspective view, for the whole domain [0, 2π] * [0, 2π]), from top to bottom the (5, 0)⊕(5, 0), (4, 1)⊕(5, 0)and (3, 2)⊕(5, 0) solutions.In Figure4can be seen the three remaining solutions (4, 1) ⊕ (4, 1), (3, 2) ⊕ (4, 1) and (3, 2) ⊕ (3, 2). All these solutions are found to be linearly stable, although all the solutions of Figure3((5, 0)⊕ something) are extremely sensitive to noise. It must be pointed out that most of these solutions would have been almost impossible to find from a time integration of the 2D Sivashinsky equation (Equation (3)) because of this sensitivity to noise, and it is likely that obtaining them from a steady version of (3) would have been very difficult too.

  , 1) ⊕ (5, 0) (3, 2) ⊕ (5, 0) (4, 1) ⊕ (4, 1), (3, 2) ⊕ (4, 1) and (3, 2) ⊕ (3, 2) solutions. Higher values of 1/ν would correspond to an increasing number of stable stationary solutions.Naturally, neither Figure1(in 1D) or Figure6(in 2D) contain all the possible stationary solutions. In 1D Guidi and Marchetti[11] have introduced the concept of interpolating solutions, which are unstable solutions connecting different branches of stable solutions in the previous bifurcation diagrams. In I, the present author has shown that this leads to a complex network of solutions, which was called web of stationary solutions. But now in two dimensions, we have the possibility, when two branches φ 1 and φ 2 exist for a parameter 1/ν to create the 2D branch φ 1 ⊕ φ 2 . This construction leads to a bifurcation diagram (with as before 1/ν < 14, i.e. not very large flames) with a truly huge number of different stationary solutions (several thousands of branches). The comparison with Figure6shows that most of these solutions are linearly unstable.The author would like to insist here on different points. First, it is only possible to obtain such an incredible number of stationary solutions because of two properties of the Sivashinsky equation, the pole decomposition, which transforms the search of stationary solutions in one dimension in a 0D problem, then the possibility to add 1D solutions in order to get 2D rectangular solutions of the Sivashinsky equation. In the Kuramoto-Sivashinsky equation case (a non linear equation with a different growth rate but the same non linear term)

Figure 1 :

 1 Figure 1: Stable stationary solutions in 1D: amplitude ∆φ vs 1/ν. All the different branches are only plotted for the values of 1/ν where they are stable. A notation like (3,2) means that 3 poles are located at x = 0, and 2 poles at x = π

Figure 2 :

 2 Figure 2: Perspective view of the (from top to bottom) (5, 0)⊕(0), (4, 1)⊕(0), (3, 2)⊕(0) stationary solutions for 1/ν = 10 and b = π/10. The solution is plotted in the interval [0, 2π] * [0, 2b] because it is easier to visualize. Actually, Neumann boundary conditions are satisfied in [0, π] * [0, b] (one fourth of the domain shown).

Figure 3 :

 3 Figure 3: Perspective view of the (from top to bottom) (5, 0) ⊕ (5, 0), (5, 0) ⊕ (4, 1), (5, 0) ⊕ (3, 2) stationary solutions for 1/ν = 10. The solution is plotted in the interval [0, 2π] * [0, 2π] because it is easier to visualize. Actually, Neumann boundary conditions are satisfied in [0, π] * [0, π] (one fourth of the domain shown).

Figure 4 :

 4 Figure 4: Perspective view of the (from top to bottom) (4, 1) ⊕ (4, 1), (4, 1) ⊕ (3, 2), (3, 2) ⊕ (3, 2) stationary solutions for 1/ν = 10. The solution is plotted in the interval [0, 2π] * [0, 2π] because it is easier to visualize. Actually, Neumann boundary conditions are satisfied in [0, π] * [0, π] (one fourth of the domain shown).

Figure 5 :Figure 6 :

 56 Figure 5: Amplitude vs time for 1/ν = 10 , starting from a (4, 1) ⊕ (3, 2) solution. A gaussian white noise (amplitude a = 0.001 ) is imposed on this solution when time is smaller than 0.5. The solution returns exponentially toward the initial solution.

Figure 7 :Figure 8 :

 78 Figure 7: Stationary solutions in 2D for a square domain : amplitude ∆φ vs 1/ν (figure with allthe solutions obtained by addition of 1D stationary solutions). When two 1D branches found in I coexist for a certain value of 1/ν, a 2D branch is created, whose amplitude is the sum of the 1D amplitudes.

Figure 9 :

 9 Figure 9: Distance to the main stationary solutions vs time for 1/ν = 10 and a = 0.01. A distance is a norm of the difference between the solution at a given time and the stationary solution. This diagram makes it possible to verify if a solution at a given time is indeed close to a stationary solution.

Figure 10 :

 10 Figure 10: Amplitude vs time for 1/ν = 10 and a = 0.1 (moderate noise amplitude). This figure suggests that the solution is often close to the (3, 2) ⊕ (3, 2) solution. It will be shown in the following figures that it is close to the (3, 2) ⊕ (4, 1) solution only for times around 50.

Figure 11 :

 11 Figure 11: Distance to the different symmetries of the (3, 2) ⊕ (3, 2) solution vs time for 1/ν = 10 and a = 0.1. The noise is sufficiently large to induce transitions between the different symmetries of this fundamental level.

Figure 12 :

 12 Figure 12: Distance to the different symmetries of the (3, 2) ⊕ (4, 1) solution (first excited level) vs time for 1/ν = 10 and a = 0.1. The solution is only close to one symmetry of the (3, 2) ⊕ (4, 1) solution for times around 50 (after apparently a transition from a (3, 2) ⊕ (3, 2) solution).

Figure 13 :

 13 Figure 13: Solution at time 120.1 for 1/ν = 10 and a = 0.1, plotted as a grayscale figure (white: minimum of φ, black: maximum of φ). The solution is plotted in the interval [0, 2π] * [0, 2π].