Surface Wave Instabilities, Period Doubling and an Approximate Universal Boundary
of Bubble Stability at the Upper Threshold of Sonoluminescence

Joachim Holzfuss
Institut für Angewandte Physik, TU Darmstadt, Schloßgartenstr. 7, 64289 Darmstadt, Germany

Phys. Rev. E 77, 066309 (2008)


Numerical results are presented for hydrodynamic instabilities surrounding the parameter space of sonoluminescing bubbles. The Rayleigh-Taylor instability is shown to limit bubble oscillations only at a small region near a border at small radii in noisy environments. Also, two different noise induced instabilities by secondary collapses are identified. A reason for period doubled, anisotropic emission is found to be due to coupling of radial and surface oscillations. Furthermore an approximate universal border is shown to exist above which the bubble is destroyed by the parametric instability. The results are compared to experimental results from several publications using different experimental setups.


1 Introduction

An air seeded bubble in water driven by a strong soundfield can oscillate stably for very long time intervals [1]. Because of enormous compression rates high temperatures occur leading to visible light emission at collapse time [2], termed sonoluminescence. A prerequisite for the stability of the process is an equilibrium set up by diffusive influx of dissolved air and outflux of dissociated reaction products. For air seeded bubbles a dynamical equilibrium is achieved due to the existence of argon as a noble gas [3]. Sonoluminescing bubbles exist only in a limited parameter range. An upper threshold is known to exist in the direction of high pressures and large ambient radii. The bubble disappears and cannot be seeded again until the driving pressure is lowered. Surface waves have been made responsible for the destruction. As the amplitude of these oscillations can increase to the size of the radial oscillation amplitude, the bubble splits into fragments [4] and may not survive this catastrophic event.

Radial and surface oscillations of bubbles are calculated and their type and behavior is determined. The total parameter space of sonoluminescing bubbles is scanned for different instabilities and the results are mapped.

2 Radial Oscillations

The Gilmore model [5] describing the radial motion of a bubble in a compressible liquid is integrated numerically.

  (      )     (       )   (      )    (      )
       R˙   3       R˙          ˙R   R       R˙
RR¨ 1- C- + 2R˙2 1 - 3C- =H  1 + C-+ C-H˙ 1- C-

    ∫ p(R,˙R)                {   }n
H =        ρ-1dp , p-+-B-=   ρl
     p∞     l      p0 + B    ρ0

           ∘ ---|        (           )
              dp||         p(R,R˙)+ B  (n-1)∕2n
C = c|r=R =   dρ-||    = c0  ---p-+-B---
               l r=R         0

p(R,R˙) =   pg(R )- 2σ-- 4ηR˙
            (      R) ( R     )γ
                 -2σ   R30 --b3
  pg(R ) =    p0 + R0   R3 - b3      (4)
R is the bubble radius, C, ρl, and p(R,) are the speed of sound in the liquid, its density, and the pressure at the bubble wall, respectively. The Tait equation is taken as the equation of state for water using n=7.025, B=3046 bars [6] as parameters. pg is the pressure in the bubble. H is the enthalpy difference of the liquid at pressure p and p(R,) at the bubble wall. p = p0 + pA cos(2πft) is the pressure at infinity, the ambient pressure p0 = 1 atm is increased by the hydrostatic pressure above the bubble. pA is the maximal ultrasound driving amplitude with frequency f = 20 kHz. The ambient bubble radius R0 is calculated at pA = 0. σ is the surface tension and η is the viscosity of the liquid. The fluid parameters are taken from tabulated values [7].

b is a van der Waals hard-core radius and γ a polytropic exponent. Its value is set between 1 (=isothermal) and the adiabatic exponent of the gas according to an instantaneous Péclet number [8910] Pe = R02|(t)|R(t)κ, reflecting thermal conduction at the involved time scales. κ is the thermal diffusivity of the gas. To avoid a change to isothermal behavior at the bubble wall turning point during maximum compression, the bubble wall velocity is kept to its maximum during positive bubble accelerations at collapse. It is smoothly approaching the real velocity during the rest of the cycle.

b and γ are updated during calculations to reflect the actual gas content. The value of the thermal diffusivity κ is variable, as it depends on the varying density ρg of the gas.

     k     ρg(R0)
κ = ρ-c-=  ρ-(R)-2⋅10-5 m2s-1
     g p    g

k is the thermal conductivity and cp the specific heat at constant pressure. The value of κ is scaled with the ratio of ambient gas density to actual density.

The density in the bubble is calculated by

ρg(R ) =   nt4otalM-                (6)
4   3     ntotalRgasTB-   4  3
3πR   =        pg     +  3πb      (7)
pg is the pressure in the bubble, Rgas = 8.3143 J mol-1K-1 is the gas constant. The van der Waals hard-core b is calculated as an average from the tabulated hard-core values and the number of moles of the i different gases, i = N2, O2, Ar, H2O. M is the molar mass averaged over all gases and ntotal the total number of moles. b and M are updated continuously.

The temperature TB is taken to be uniform within the bubble. It is calculated via the adiabatic compression of a van der Waals gas by

       (   3   3)γ-1
TB = T0  R-03 --b3
         R  - b

with the ambient liquid temperature T0. Other more complex approaches exist that include thermal conduction and a temperature jump across the bubble interface [1112].

Evaporation and condensation of water molecules at the bubble wall [121331415] are included in the model for the bubble dynamics, as experimental results [1617] stress the importance of a decrease of the polytropic exponent induced by water vapor at bubble collapse. A simple Hertz-Knudsen model [1819] for the change of moles nH2O of water vapor in the bubble is

n˙H2O   =  ˙neHv2aOp- ˙ncHon2dO
           4πR2 α¯c(T )(              )
       =  ----------s- ρsagt,H2O - Γ ρg,H2O  .(9)
          MH2O    4
α is the constant evaporation coefficient (also called accommodation coefficient, sticking probability),
      ∘ --------
¯c(T ) =   8RgasTs-
  s      πMH2O

is the average velocity of molecules of a Maxwell-Boltzmann distribution, ρg,H2O is the density of water vapor of molar weight MH2O in the bubble, ρg,H2Osat the saturated vapor density [7]. The bubble surface temperature is taken as Ts = T0. The density of water vapor ρg,H2O depends on the bubble dynamics and is calculated along with the bubble equation with the help of (6). The simple model (9) takes the temperature distributions in the bubble and liquid as fixed and does not capture all effects occurring during evaporation, as would do more complex treatments [202122]. Γ is a correction factor for nonequilibrium conditions induced by mass motion of vapor and bubble wall movement (Schrage correction) [231114].

               (       ∫        )
     -Ω2    √--     -2-  Ω - x2
Γ = e   - Ω  π  1 - √π- 0 e   dx

    R˙- v
Ω = -----H2O

Ω is a ratio of velocities, is the bubble wall velocity, vH2O is the vapor velocity and cpeak is the velocity belonging to the peak of the Maxwell-Boltzmann velocity distribution.

          ∘ --------
c   (T) =   2RgasTs
peak  s      MH2O

The change of mass per unit time and unit area j can be expressed as

j = n˙H2OMH2O- = ρg,H2O ( ˙R - vH2O )

and when inserted into (12) leads to

Ω =  nH2O 3cpeak

for a spherical bubble volume. For small values of Ω the approximation Γ = 1 - Ω√ π- can be made [23] leading to an effective constant evaporation coefficient of αeff = 2α∕(2 - α) in (9) and an effective correction factor Γeff of unity. Calculations show, that Γ varies by as much as 20% around unity during the collapses. However, in the observed parameter range almost no notable difference exists in the amount of water vapor at collapse time between results of calculations using a constant value of αeff and setting Γeff = 1, and of those using (9) with the variable expressions (11) and (15) together with a smaller value of the evaporation coefficient α = 2αeff(2 + αeff). Therefore a constant value of αeff = 0.4 [24] is taken in the following calculations.

3 Surface Waves

A bubble in a liquid is not only oscillating radially, but the surface may also exhibit a wavy oscillation as has been observed earlier [25]. The spherical stability of a bubble has been analyzed by a great number of authors [2628292730]. The problem is involved since during one oscillation the direction of bubble acceleration changes as well as the density within the bubble.

A linear analysis is made studying the perturbation of the bubble shape [3126]

r(t,Θ,φ) = R(t) + an(t)Ymn (Θ, φ) ,

where R(t) is the instantaneous bubble mean radius, Y nm a surface harmonic of degree n and order m and an its amplitude. Prosperetti developed a theory for surface oscillation of a bubble [32] leading to an integro-differential equation for the amplitude disturbance an. In the linear analysis, the perturbation dynamics is independent of m. For small liquid viscosity integral terms can be dropped. Neglecting density and viscosity of the gas content in a liquid with constant density ρl = ρ0 a simpler expression is [33]

       [                       ]
         R˙                --η-
¨an  +   3R  +2(n + 2)(2n +1)ρlR2  ˙an
    +  (n- 1)  - ¨R-+ (n + 1)(n + 2)-σ--
                 R               ρlR3
                      +2 (n+ 2)-η ˙R an = 0(17)
For a bubble of fixed radius R, the time derivatives vanish and the natural frequency wn and the damping coefficient bn for the n-th mode can easily be deduced from (17) as
 2                       --σ-
ωn  =  (n - 1)(n +1)(n+ 2)ρlR3    (18)
βn  =  (n + 2)(2n +1)ρlR2         (19)
The first natural frequencies ωn2π, n = 2,3,4 can be calculated with the usual parameters to 2.638, 4.817 and 7.226 MHz for a bubble of 5 μm ambient radius and 3.687, 6.731 and 10.097 MHz for a 4 μm bubble. This suggests that the n=2-surface mode is the first excited mode due to its lower natural frequency.

Lohse & al. [3435] went one step beyond the above approximation and arrived at a boundary layer type approximation, where the vorticity only acts within a small distance from the bubble. Using their BLA approximation on the original equations [32] including the dependence on the time varying gas density ρg [32] and neglecting gas viscosity [3637] one gets

 ¨an  +  Bn ˙an + Anan = 0         (                               )                          (20)
         R˙  (  ρl    ρg) -1  2η   n(n + 2)2   1
Bn   =  3R-+   n-+1-+ -n    × R2-  -n-+-1--1+-2δ∕R-- (n- 1)(n+ 2)                           (21)
        (          )
          -ρl--  ρg  -1
An   =    n+ 1 + n                                                                          (22)
          [(                ) ¨                       ˙ (                                   ) ]
        ×    n+-2ρg - n--1ρl  R-+ (n - 1)(n + 2) σ-+2η R-- (n- 1)(n+ 2)- n(n---1)(n-+-2)---1----
              n       n+ 1    R               R3     R3                     n+ 1     1+ 2δ∕R
where δ is the boundary layer thickness, suggested [38] to be defined as
       (∘ -η-- R )
δ = min    ---,---
           ρlω  2n

where ω∕2π is the ultra sound driving frequency. In the following calculations the boundary layer is set to a value δ = 0, which, as can be seen by the results, models stability boundaries in a fairly good agreement with experimentally measured boundaries [37].

3.1 Rayleigh-Taylor Instability

Basically two types of instability have been identified. The Rayleigh-Taylor instability (RT) [3940] acts on very short time scales in the order of the duration of the bubble collapse. It only acts when the bubble wall acceleration ¨R is positive such that the lighter gas inside the bubble is accelerated into the heavier fluid, namely water. Small disturbances an R may then be amplified exponentially, eventually leading to a bubble breakoff. This is most prominently taking place around the main bubble collapse. Arguments exist in the current literature that the RT instability may be suppressed by a proper account of thermal effects in the gas [2833]. A suppression when comparing to results in [34] is also calculated [37], when the increasing density of the bubble’s interior during the collapse is accounted for (eqs. 20-23).


Figure 1: Rayleigh-Taylor instability: The driving pressure is 1.55 bars, the mean ambient bubble radius 0.903 μm and the driving frequency is 20 kHz. Shown is the ratio of amplitude of surface distortion to mean radius. Maxima and minima are seen around main bubble collapse and subsequent afterbounces. The n = 2 (upper graph) surface harmonic shows a smaller distortion at the main collapse which is amplified by the afterbounces. The n = 6 (lower graph) surface harmonic shows a large distortion at the main collapse and hardly any reaction to the afterbounces.

Different methods for the calculation of the RT instability have been developed. Sudden growth of the surface oscillation amplitude to the size of the bubble radius defines the instability threshold resulting in a breaking of the bubble. To induce this instability numerically, different strategies exist. In [34] a random displacement of size 0.1 nm is added as microscopic fluctuations to an in (20) during integration. Others have added a fixed initial kick to the surface oscillation amplitude of 10 nm shortly before the bubble collapse [2937], added gaussian distributed noise with 0.1 nm standard deviation as molecular fluctuations [41] , or reformulated the problem as a stochastic differential equation [36]. Here, gaussian distributed noise with 0.5 nm standard deviation is added to (20) during integration. An important quantity is the absolute value of an(t)∕R(t): A maximum exceeding a value of unity implies a splitting of the bubble. Calculation of the above ratio shows that it depends on the amount and type of noise. Changing the strength of the fluctuation by an order of magnitude is reflected as such and hinders the determination of the exact position of the RT instability in parameter space. As the amplification rate of disturbances induced by the RT instability is inversely proportional to their wavelength, higher order surface harmonics are used for the calculation. As a side effect, the afterbounce instability (see next paragraph) as a parametric effect is disappearing for high order modes (Fig. 1). The Gilmore equations together with (20) - (23) are used in Fig 1 for an argon-water vapor bubble with van der Waals hard-core driven at 20 kHz including water evaporation and adiabatic-isothermal switching according to the Péclet number. A drastic increase of the amplitude of the n = 6 mode to a value larger than the threshold of breakup is seen during main bubble collapse.


Figure 2: (Color online) Rayleigh-Taylor instability of an argon-water vapor bubble: In the parameter space spanned by driving pressure and equilibrium radius the average maximum of the ratio of instantaneous surface deformation to mean radius is color coded. The maximum deformation ratio of modes n = 2 to 6 is taken around the main collapse and averaged over ten consecutive driving oscillation cycles. Within the contour lines the value is larger than unity implicating splitting of the bubble. The equilibrium radius is taken as the value of radius with no periodic forcing amplitude applied.


Figure 3: Two types of afterbounce instabilities: Upper: Driving pressure is 1.4 bars, ambient radius 5.23 μm. Lower: Driving pressure is 1.7 bars, ambient radius 2.51μm. Driving frequency is 20 kHz. The times of the main collapses are indicated by arrows.


Figure 4: (Color online) Amplification of surface mode oscillations of an argon-water vapor bubble by afterbounces: In the parameter space spanned by driving pressure and equilibrium radius the average maximum of the ratio of surface deformation (n = 2) amplitude to mean radius is color coded. The average maximum occurring during a whole driving oscillation period of 6 different cycles is taken. The individual cycles are non consecutive and randomly disturbed. Within the contour lines the value is larger than unity implicating splitting of the bubble. Two main regions of this instability can be seen to exist.

Parameter space calculations show that the existence of the RT instability is limited to a small region in phase space (Fig. 2). The radius for an undriven bubble is taken as the equilibrium radius. Due to usage of the Gilmore model that includes a better modeling for acoustic damping when compared to simpler models, the inclusion of the density correction in (20) - (23) and water evaporation-condensation it is seen, that the RT-unstable region is very small and confined to high pressures and small mean equilibrium bubble radii. The region gets smaller at increasing pressures (compare [34373642]). Furthermore, the amount of added noise in the calculation has to be increased by a factor of 5 to show a notable effect. It may be asked if the added amount of 0.5 nm average gaussian noise amplitude is too large for experimental conditions. If the answer is yes then a destruction of the bubble by a genuine Rayleigh-Taylor instability is not happening.

Figure 5: Parametric instability of an argon-water vapor bubble. The maximal amplitude of surface deformation increases during several periods of the driving sound to a value larger than the radius thus leading to bubble breakup. Left: The upper trace shows bubble radius, while the lower trace shows surface oscillation amplitude (n = 2). Equilibrium radius is R0=5.74 μm, 1.4 bars is the sound pressure. Right: Excerpt of left graph around point of splitoff, f = 20 kHz is the driving frequency for all calculations.

3.2 Instability by Secondary Collapses

An amplification of the surface mode oscillations by repeated secondary collapses (afterbounces) subsequent to the main collapse has been termed ”afterbounce instability” [34]. It may induce very large surface oscillations leading to the destruction of a bubble within a single cycle. The instability is most predominant for low order surface harmonics. Two types of this instability can be observed (Fig. 3): It can be triggered by an initial RT instability (Fig. 3, lower graph) which induces a large surface deformation not leading to breakup during the main collapse. The deformation is not damped out until the next afterbounce thus leading to an even larger surface deformation during the first afterbounce collapse while having the same short time behavior as an RT instability. This may repeat until the afterbounces finally get too small for further amplification. This type of instability can only be demonstrated numerically with a large amount of noise fluctuations of 0.5 nm standard deviation of gaussian distributed noise.


Figure 6: (Color online) Parametric stability of a argon-water vapor bubble: In the parameter space spanned by driving pressure and equilibrium radius the growth rate of surface deformation per period (n = 2) is color coded. Within the white contour lines the value is larger than unity implicating parametric growth of the surface deformation.

Another type of amplification by afterbounces is seen to start to exist near the bifurcation to parametric instability (Fig 3, upper graph): Due to the added noise the surface oscillations get very large and may exceed the minimum of the mean bubble radius thus leading to destruction. This happens only at later afterbounces which resonantly amplify the surface deformations. These deformations only grow from cycle to cycle once the parameters are set to values above the threshold for parametric instability. This type of afterbounce instability also exists for small amounts of added noise.


Figure 7: Parametric growth of surface deformation of an argon-water vapor bubble: The time dependence of the (n = 2) surface oscillation per radius ratio is shown for 4 consecutive time sections around the main collapse (marked by arrows). Data from two different parameter settings have been taken. First row: equilibrium radius R0=5.839 μm at 1.35 bars, second row: R0=6.0 μm at 1.35 bars. Note the alternating symmetry in the first row.

3.3 Parametric Instability

The parametric instability, in planar geometry also known as the appearance of Faraday waves [4344], describes the growth of the absolute magnitude of the surface deformation amplitude an over many cycles of the driving sound. If small disturbances an R are vanishing from period to period, the bubble is said to be parametrically stable. Some debate has occurred concerning the role of viscosity concerning the parametric instability [4546] resulting in an underestimation of the excitation threshold of bubble radii for low driving pressures. A semianalytical approach [29] and more exact numerical models [4748] finally showed the validity of eqs. 20 - 23 at least in the parameter region for sonoluminescence.

The parametric instability is calculated in principle using Floquet’s theorem by determing the amplification rate of a small perturbation after one period of a steady-state radial oscillation [34]. In the parameter region analyzed here this can be relaxed to monitoring the growth of the surface deformation amplitude an over several cycles. It is found that the n = 2 mode has the lowest threshold for a per-period-amplification-rate larger than unity. Figure 5 shows a typical example of parametric amplification of an n = 2 mode. During a single period the time trace shows the development of a surface oscillation. It is seen that the surface oscillations are driven by the afterbounces. The oscillation period equals one half of the driving afterbounces. The envelope shows an increase of the oscillation amplitude followed by a subsequent decrease. The remaining amplitude of the surface deformation at the time of the next main collapse is amplified by the RT instability again resulting in a parametric oscillation. If the maximal amplitude is increasing over the timescale of several periods this leads to bubble breakup.

Figure 6 shows the value of the growth rate per driving period of the n=2 surface oscillation mode in the parameter space of SBSL. Parameters of equilibrium radius and driving pressure from within or above the white contour lines show parametric instability leading to bubble breakup. As a tendency less driving amplitude and smaller bubbles are more stable. However a fine structure is superimposed onto the graph showing individual cells separated by horizontal and vertical lines. At the lines the amplification rates are smaller compared to those within each cell. This is true below the stability threshold as well as above it. Near the boundary a more complicated scenario exists: Eventually a cell contains an island surrounded by a white contour line denoting the presence of an unstable region within an otherwise stable cell.

Figure 7 shows simulations of the amplitude of the n = 2 surface oscillation mode with two parameter settings from within adjacent cells of Fig. 6. Time sections of blowups of subsequent oscillations are shown. Both rows show a parametric amplification of the surface oscillation amplitude. However, in contrast to the second, the first row shows an alternating sign symmetry of the waveform of consecutive periods. This also results in an opposite sign of the amplitude during the main collapse (marked by arrows). Whenever a parameter setting is picked from an adjacent cell, a transition from non-alternating to alternating behavior is made. This is true within the whole parameter space, for parametrically stable or unstable settings. An alternating parametrically amplified waveform is also seen in Fig. 5.

The amplitude of the n = 2 oscillation is about 2% of the mean radius at main collapse at the time the ratio of n = 2 oscillation to mean radius exceeds 1 during the afterbounce oscillations implicating breakup of the bubble.

4 Discussion

The region of hydrodynamical stability of driven bubbles in the parameter range of sonoluminescence has boundaries in the direction of increasing equilibrium radius and increasing driving pressure. The Rayleigh-Taylor instability is numerically shown to be limited to a region of very small bubbles at high driving pressures. The instability set up by secondary collapses limits the accessible parameter space also in the direction of larger equilibrium radii as well as higher pressures. The determination of the exact position of the boundaries are hindered by the amount of added noise used in the numerical procedure. When comparing to experimentally determined boundaries, the experimentally existing noise should also be considered.

The parametric instability limits the parameter space mostly in the direction of growing radii. A consequence for a sonoluminescing bubble driven near the threshold of parametric stability is that an appreciable distortion (e.g. prolate-oblate ellipsoid) exists during the collapse. This may result in a directional radiation of sonoluminescence light output [49]. This has been seen previously in experiments [505152]. Measurements presented in [515352] have shown that during unstable SBSL, where the equilibrium radius of the bubble and also its spatial position and driving pressure changes, a bubble may show different radiation behavior, from uniform to non-uniform and vice-versa. This may be understood as radiation from a bubble entering and leaving islands of symmetrically different parametric stability. Interestingly, measurements in [5152] show a double periodic cycle in the intensity of spatially directed radiation near the upper threshold of SBSL and during unstable SBSL. This can be well understood with the results of Fig. 7, where the surface amplitude at collapse alternates for subsequent periods of the driving revealing a symmetry change of the bubble surface at main collapse with a doubled period. As the average radial oscillation of a bubble is only driving the surface oscillation and experiences no back coupling, it would be unaffected in this scenario, resulting in a constant periodic collapse as measured in [52]. The absolute value of the parametric stability is only slightly above one in the first 2 (for increasing radius) elongated cells above the threshold in Fig. 6. Thus it may be possible that nonlinear surface mode effects may inhibit the growth leading to bubble destruction. Also in the case of a bubble growing by diffusion [4] it may pass through an unstable island fast enough without getting destructed [5352].

The existence of the cells in Fig. 6 is a consequence of coupling of nonlinear oscillators (one sinusoidally driven radial oscillator driving the surface oscillator) and does not depend on the types of models and approximations used. Signs of a cell-like structure can also be seen in graphs showing parametric boundaries published by other authors using different models and approximations [383327]. The exact position of lines of instability in parameter space spanned by ambient radius and driving pressure do depend somewhat on the models and parameters used (in accordance with the current literature). Their selection can only be justified by comparison with experimental data. Models using a nonuniform bubble interior should find inner high pressure peaks or shock waves being detached from the elliptical bubble surface resulting in a nonspherical core. This will lead to a non isotropic emission as well. In fact, molecular dynamics simulations have been published [54] that report non isotropic emission upon collapse of slightly ellisoid bubbles during SBSL condistions.

An instability that has as yet not been analyzed in SBSL research but may be of importance, is the Richtmyer-Meshkov instability [555657]. It occurs when a shock wave interacts with the gas-fluid interface. The Kelvin-Helmholtz instability [57] occurs when there is a fluid motion parallel to the interface between gas and water, which may be present in the form of microstreaming [5859].


Figure 8: Maximal bubble radius times driving frequency at parametric instability threshold as a function of ratio of driving pressure over ambient pressure. Bubble oscillations below the line are parametrically stable. Calculations are shown for ambient pressures of +/-1.4 kPa (square,triangle) around 1 atm (circle) with a driving frequency of 20 kHz (filled symbols), 14 kHz (hollow square), 24 kHz (star) and 30 kHz (hollow circle) at 1 atm. The boundary layer is δ = 0. The line shows a linear regression over all data points.


Figure 9: Maximal bubble radius times driving frequency as a function of normalized pressure of various data obtained from literature: cross [6061], star [62], filled circles [63], down triangles [16], up triangle [64], filled squares: own data, hollow circle: [65]. Where indicated the ambient pressure (set to 1 atm if unknown) has been augmented by hydrostatic pressure guessed from published cell dimensions. The line is determined from Fig. 8.

5 Approximate Universal Boundary

The calculated threshold for the parametric instability using the boundary layer approximation in (23) is found to be lower than the threshold at which a bubble disappears in experiments. It has been shown [37] that smaller values of the size of the boundary layer δ (eq. 23) result in a shift of the parametric instability to larger ambient radii. Actually, setting δ = 0 resulted in the best coincidence of numerical and experimental stability boundaries. One has to cope with the fact that it is somewhat difficult to measure the bubble’s ambient radius. This is because of limited optical resolution and the question, when and how to measure an ambient radius in the case of sinusoidal driving. A solution is to measure the maximal radius Rmax of the bubble during an oscillation.

Fig. 8 shows RmaxPIf, which is the maximal bubble radius at the first lower border of parametric instability times the driving frequency, as a function of the driving pressure over ambient pressure ratio. An almost linear dependence is suggested over a large parameter variation covering the whole SBSL driving parameter space. A linear regression over all data points results in RmaxPIf = -2.296 + 2.648(pA∕p0) with a standard error of 0.01 in both constants. A more detailed calculation (Fig. 11) of the boundary shows the same but very squeezed roughness as Fig. 6, with all the long island lines running parallel to the line in Fig. 8.


Figure 10: Maximal bubble radius times driving frequency at parametric instability threshold as a function of normalized pressure. Data from various models and approximations are shown. Circles are calculated from the Gilmore bubble model using all equations specified earlier. Crosses denote the Keller-Miksis model (KM). The modified Rayleigh-Plesset bubble model (RPm) is shown with several changes applied to demonstrate the robustness of the approximation of a universal boundary. The line is determined from Fig. 8.

Various data of the maximal bubble radius at the border of SBSL extinction have been obtained from the literature and are shown in Fig. 9. The data sets have been taken by different authors, with different setups, driving frequencies and pressures. Ambient pressures have not always been published, nor has the hydrostatic pressure induced by the water column on top of the bubble, which has hence been set to an elaborated estimate. A very good agreement with the calculated linear relationship in Fig. 9 is seen. In particular the detailed results from [63] agree well. The outmost lying circle corresponds to a published ambient pressure being somewhat ”out of sequence”. The down triangles (from results of Fig. 5, curves a, b in [16]) lie outside because they are taken at an ambient pressure of 2 bars, which is not in the applicability range of the dependence suggested by Fig. 8. The results represented by the squares stem from a rather noisy experiment. The hollow circle is a maximally driven nitrogen bubble [6665] at 9 deg C close to breakup.

Experimental and numerical studies have been made concerning the dependence of light output and change of ambient pressure [6768]. The authors find a linear increase of the number of photons as the ambient pressure p0 is decreased. Here we are concerned about the maximally obtainable bubble radius close to the border of instability. At constant frequency and driving pressure a decreasing of p0 increases the maximal bubble radius at the border of parametric instability too, eventually leading to a higher number of emitted photons.

The nearly perfect linear dependence of the threshold line RmaxPIf as a function of normalized pressure pA∕p0 in the parameter range for SBSL suggested by Figs. 8 and 9 can be motivated: A larger maximal radius leads to a more violent bubble collapse. The resulting larger amplitudes of afterbounces lead to stronger excitation of surface waves. An increasing frequency of the collapses leads to a further amplification as the surface waves have less time to be damped between each excitation. The dependent axis has units of a velocity. Parametric instability is induced when this threshold velocity is passed over. The requirement of the pressure axis normalization by the ambient pressure is also seen in the case of ambient pressure dependence of diffusive stability lines [69703]. A larger ambient pressure decreases the maximum radius during the oscillation.

The existence of an approximate universal boundary of the parametric instability is shown to be robust against changes of the bubble equation and various approximations made. Fig. 10 shows parametric instability lines for the Gilmore, Keller-Miksis [71] and modified Rayleigh-Plesset equation [3] driven at f = 20 kHz. The equations are augmented with the other model equations mentioned earlier. Hardly any difference on the choice of bubble equation is visible. Due to its more complete modeling capabilities of radiated sound, the Gilmore equation shows slightly smaller values of the maximal radius for large forcing amplitudes due to acoustic energy radiation upon collapse in the form of shock waves. Calculations for the modified Rayleigh-Plesset model with an additional change in approximation have been made: A purely adiabatic bubble compression with constant adiabatic exponent (eq. (8)) shows the largest deviation from the line calculated from Fig. 8 to values that are also well below the experimental data points from Fig. 9. The neglection of the evaporation-condensation equations also results in overall smaller values, as does the inclusion of the BLA approximation via eq. 23. The latter is in accordance with the results from [37] on the ambient radius. The results for a different equation of state for the gas inside the bubble (ideal gas) do not show any considerable difference to the results for a van der Waals hard-core gas. This motivates the assumption that using a more realistic, but computationally expensive soft shell model [72] also leads to a very small difference between the results.

Simulations using molecular dynamics have been published [547374] that use a bubble equation as a launch condition. Backreactions of non isotropic bubble interior to surface oscillations have not yet been studied in this context. Analysis of those may lead to further agreement of theory and experiment.

Fig. 11 shows the regions of all instabilities analyzed before in a single graph. In comparison it is seen that the threshold for the noise induced after bounce instability AB1 is shifted to smaller maximal radii below the parametric instability. The Rayleigh-Taylor instability RT occurs at smaller maximum radii near the boarder of the afterbounce instability AB2. A small gap with none of the above instabilities extends to larger maximal radii with increased driving pressure.


Figure 11: (Color online) Phase space spanned by driving pressure and maximal radius showing regions of Rayleigh-Taylor (RT)-, afterbounce (AB)- and parametric (PM) instabilities. A region for high forcing amplitudes exists between AB1 and AB2 where none of the above instabilities exist.

6 conclusions

Nonisotropic radiation in combination with spatially directed period doubling has been observed in experiments by M. Levinsen & al.. In order to find an explanation of this, different hydrodynamical instabilities have been mapped throughout the whole parameter space. The parametric instability induces slowly growing surface waves, calculated as an prolate-oblate oscillation. The structure of the growth rate of this instability mapped in parameter space shows a ragged threshold line, consisting of multiple islands being separated by stability zones. The dynamical behavior between adjacent islands is different. In one islands the surface oscillations catch up on each other such that the bubble is always prolate or always oblate at collapse. In the neighboring island a period 2 oscillation is seen such that at one period the bubble is oblate while in the next it is prolate. This last type of oscillation will result in a spatially directed emission of light with a period doubled intensity just as observed in the experiment.

The Rayleigh-Taylor instability occurring around collapse upon acceleration of the lighter fluid (gas) of the bubble into the heavier fluid (water) is calculated using varying gas density and higher order disturbances. This instability is found to occur only at small ambient radii and stronger forcing and should therefor not occur at the upper driving pressure threshold of SBSL. Two different types of an afterbounce instability have been identified and shown to exist in different regions of the parameter space. For experiments with a high amount of noise (sound disturbances, higher ambient temperature) one type of afterbounce instability may impose an upper limit to SBSL.

A global threshold function has been found in the driving pressure - maximum radius space at the parametric instability. Within the limits of the ragged islands mentioned earlier, which are visible only in very detailed calculations, a rescaling of the pressure axis with the ambient pressure and the maximum radius at threshold with the driving frequency reveals a straight line dependence in the SBSL parameter range. All available experimental data from literature have been inserted into a graph. Within experimental accuracy a global threshold value of the parametric instability as a function of rescaled pressure can be deduced.

The author thanks M. Levinsen.


[1]   D. F. Gaitan, L. A. Crum, C. C. Church, and R. A. Roy, J. Acoust. Soc. Am. 91, 3166 (1992).

[2]   B. P. Barber, R. A. Hiller, R. Löfstedt, S. J. Putterman, and K. R. Weninger, Phys. Rep. 281, 65 (1997).

[3]   M. P. Brenner, S. Hilgenfeldt, and D. Lohse, Rev. Mod. Phys. 74, 425 (2002).

[4]   J. Holzfuss, Phys. Rev. E 71, 026304 (2005).

[5]   F. R. Gilmore, Report 26-4, The growth or collapse of a spherical bubble in a viscous compressible liquid, Hydrodynamics Laboratory, California Institute of Technology, Pasadena, Cal., USA (unpublished), (1952).

[6]   J. Holzfuss, M. Rüggeberg, and A. Billo, Phys. Rev. Lett. 81, 5434 (1998).

[7]   At 20 deg C the following parameters have been taken in the calculations: surface tension σ = 0.0725 N m-1, viscosity η = 0.001 N s m-2, saturated vapor density ρg,H2Osat = 0.0173 kg m-3, diffusion constant DAr = 2.0 × 10-9 m2s-1, diffusion constant DO2 = 1.9 × 10-9 m2s-1, diffusion constant DN2 = 1.9 × 10-9 m2s-1, solubility Ar = 0.061 kg m-3, solubility O2 = 0.0444 kg m-3, solubility N2 = 0.020 kg m-3, water density ρ0 = 998.23 kg m-3.

[8]   D. Lohse and S. Hilgenfeldt, J. Chem. Phys. 107, 6986 (1997).

[9]   A. Prosperetti, J. Acoust. Soc. Am. 61, 17 (1977).

[10]   M. S. Plesset and A. Prosperetti, Ann. Rev. Fluid Mech. 9, 145 (1977).

[11]   K. Yasui, J. Acoust. Soc. Am. 98, 2772 (1995).

[12]   K. Yasui, J. Phys. Soc. Jap. 66, 2911 (1997).

[13]   V. Kamath, A. Prosperetti, and F. N. Egolfopoulos, J. Acoust. Soc. Am. 94, 248 (1993).

[14]   B. D. Storey and A. J. Szeri, Proc. R. Soc. London A 456, 1685 (2000).

[15]   B. D. Storey and A. J. Szeri, Phys. Rev. Lett. 88, 074301 (2002).

[16]   G. E. Vazquez and S. J. Putterman, Phys. Rev. Lett. 85, 3037 (2000).

[17]   R. Toegel, B. Gompf, R. Pecha, and D. Lohse, Phys. Rev. Lett. 85, 3165 (2000).

[18]   H. Hertz, Annalen der Physik 17, 178 (1882).

[19]   M. Knudsen, Kinetic theory of gases (Methuen and Co., London, 1950).

[20]   M. Bond and H. Struchtrup, Phys. Rev. E 70, 061605 (2004).

[21]   D. J. E. Harvie and D. F. Fletcher, J. Heat Transfer 123, 486 (2001).

[22]   G. Hauke, D. Fuster, and C. Dopazo, Phys. Rev. E 75, 066310 (2007).

[23]   R. W. Schrage, A theoretical study of interphase mass transfer (Columbia University Press, New York, 1953).

[24]   I. W. Eames, N. J. Marr, and H. Sabir, Int. J. Heat Mass Transfer 40, 2963 (1997).

[25]   M. Kornfeld and L. Suvorov, J. Appl. Phys. 15, 495 (1944).

[26]   H. W. Strube, Acustica 25, 289 (1971).

[27]   Y. An, T. Lu, and B. Yang, Phys. Rev. E 71, 026310 (2005).

[28]   A. Prosperetti and Y. Hao, Phil. Trans. R. Soc. Lond. A 357, 203 (1999).

[29]   V. A. Bogoyavlenskiy, Phys. Rev. E 62, 2158 (2000).

[30]   A. A. Aganin and T. S. Guseva, J. Appl. Mech. and Techn. Phys. 46, 471 (2005).

[31]   M. S. Plesset, J. Appl. Phys. 25, 96 (1954).

[32]   A. Prosperetti, Q. Appl. Math. 34, 339 (1977).

[33]   Y. Hao and A. Prosperetti, Phys. Fluids 11, 1309 (1999).

[34]   S. Hilgenfeldt, D. Lohse, and M. P. Brenner, Phys. Fluids 8, 2808 (1996).

[35]   S. Hilgenfeldt, D. Lohse, and M. P. Brenner, Phys. Fluids 9, 2462 (1997).

[36]   U. H. Augsdörfer, A. K. Evans, and D. P. Oxley, Phys. Rev. E 61, 5278 (2000).

[37]   L. Yuan, C. Y. Ho, M.-C. Chu, and P. T. Leung, Phys. Rev. E 64, 016317 (2001).

[38]   M. P. Brenner, D. Lohse, and T. F. Dupont, Phys. Rev. Lett. 75, 954 (1995).

[39]   G. I. Taylor, Proc. R. Soc. London, Ser. A 201, 192 (1950).

[40]   D. J. Lewis, Proc. R. Soc. London, Ser. A 202, 81 (1950).

[41]   H. Lin, B. D. Storey, and A. J. Szeri, Phys. Fluids 14, 2925 (2002).

[42]   H. Lin, B. D. Storey, and A. J. Szeri, J. Fluid Mech. 452, 145 (2002).

[43]   M. Faraday, Philos. Trans. R. Soc. London 121, 299 (1831).

[44]   L. Rayleigh, Phil. Mag. 16, 50 (1883).

[45]   S. J. Putterman and P. H. Roberts, Phys. Rev. Lett. 80, 3666 (1998).

[46]   M. P. Brenner, T. F. Dupont, S. Hilgenfeldt, and D. Lohse, Phys. Rev. Lett. 80, 3668 (1998).

[47]   P. H. Roberts and C. C. Wu, Phys. Fluids 10, 3227 (1998).

[48]   C. C. Wu and P. H. Roberts, Phys. Lett. A 250, 131 (1998).

[49]   A. Madrazo, N. García, and M. Nieto-Vesperinas, Phys. Rev. Lett. 80, 4590 (1998).

[50]   K. Weninger, S. J. Putterman, and B. P. Barber, Phys. Rev. E 54, R2205 (1996).

[51]   M. T. Levinsen, N. Weppenaar, J. S. Dam, G. Simon, and M. Skogstad, Phys. Rev. E 68, 035303(R) (2003).

[52]   J. S. Dam and M. T. Levinsen, Phys. Rev. Lett. 94, 174301 (2005).

[53]   J. S. Dam, M. Levinsen, and M. Skogstad, Phys. Rev. E 67, 026303 (2003).

[54]   B. Metten and W. Lauterborn, in Nonlinear Acoustics at the Turn of the Millennium, ISNA 15, edited by W. Lauterborn and T. Kurz (AIP, New York, 2000), pp. 429–432.

[55]   R. D. Richtmyer, Commun. Pure Appl. Math. 13, 297 (1960).

[56]   E. E. Meshkov, Fluid Dyn. 4, 101 (1969).

[57]   S. Chandrasekhar, Hydrodynamic and hydromagnetic stability (Dover Publications Inc., New York, 1981).

[58]   T. G. Leighton, The Acoustic Bubble (Academic Press, London, 1994).

[59]   T. Verraes, F. Lepoint-Mullie, T. Lepoint, and M. S. Longuet-Higgins, J. Acoust. Soc. Am. 108, 117 (2000).

[60]   D. Krefting, R. Mettin, and W. Lauterborn, Phys. Rev. Lett. 91, 174301 (2003).

[61]   R. Mettin, the upper limit of the driving pressure is given by 1.34 bars. Private communication.

[62]   D. F. Gaitan and R. G. Holt, Phys. Rev. E 59, 5495 (1999).

[63]   C. R. Thomas, R. A. Roy, and R. G. Holt, Phys. Rev. E 70, 066301 (2004).

[64]   B. P. Barber, C. C. Wu, R. Löfstedt, P. H. Roberts, and S. J. Putterman, Phys. Rev. Lett. 72, 1380 (1994).

[65]   J. Holzfuss and M. Levinsen, Phys. Rev. E 77, 046304 (2008).

[66]   M. T. Levinsen and J. S. Dam, Europhys. Lett. 80, 27004 (2007).

[67]   L. Kondic, C. Yuan, and C. K. Chan, Phys. Rev. E 57, R32 (1998).

[68]   M. Dan, J. D. N. Cheeke, and L. Kondic, Phys. Rev. Lett. 83, 1870 (1999).

[69]   A. Eller and H. G. Flynn, J. Acoust. Soc. Am. 37, 493 (1965).

[70]   R. Löfstedt, K. R. Weninger, S. Putterman, and B. P. Barber, Phys. Rev. E 51, 4400 (1995).

[71]   J. B. Keller and M. Miksis, J. Acoust. Soc. Am. 68, 628 (1980).

[72]   G. A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Oxford University Press, New York, 1998).

[73]   S. J. Ruuth, S. Putterman, and B. Merriman, Phys. Rev. E 66, 036310 (2002).

[74]   A. Bass, S. Putterman, B. Merriman, and S. J. Ruuth, J. Comp. Phys. 227, 2118 (2008).