PHYSICAL REVIEW E 77, 066309 共2008兲

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ßgartenstraße 7, 64289 Darmstadt, Germany
共Received 31 October 2007; revised manuscript received 18 April 2008; published 17 June 2008兲
Numerical results are presented for hydrodynamic instabilities surrounding the parameter space of sonolu-
minescing bubbles. The Rayleigh-Taylor instability is shown to limit bubble oscillations only in a small region
near a border at small radii in noisy environments. Also, two different noise-induced instabilities by secondary
collapses are identified. The reason for period-doubled, anisotropic emission is found to be the 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.

DOI: 10.1103/PhysRevE.77.066309 PACS number共s兲: 78.60.Mq

I. INTRODUCTION 2␴ 4␩
p共R,Ṙ兲 = pg共R兲 − − Ṙ,
An air-seeded bubble in water driven by a strong sound R R
field 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 共SBSL兲. A prerequisite for the sta-
冉
pg共R兲 = p0 +
2␴
R0
冊冉 R30 − b3
R3 − b3
冊 ␥
. 共3兲

bility of the process is an equilibrium set up by diffusive Here, R is the bubble radius, C, ␳l, and p共R , Ṙ兲 are the speed
influx of dissolved air and outflux of dissociated reaction of sound in the liquid, its density, and the pressure at the
products. For air-seeded bubbles a dynamical equilibrium is bubble wall, respectively. The Tait equation is taken as the
achieved due to the existence of argon as a noble gas 关3兴. equation of state for water, using n = 7.025 and B
Sonoluminescing bubbles exist only in a limited parameter = 3046 bars 关6兴 as parameters. pg is the pressure in the
range. An upper threshold is known to exist in the direction bubble. H is the enthalpy difference of the liquid at pressure
of high pressures and large ambient radii. The bubble disap- p⬁ and p共R , Ṙ兲 at the bubble wall. p⬁ = p0 + pA cos共2␲ ft兲 is
pears and cannot be seeded again until the driving pressure is the pressure at infinity; the ambient pressure p0 = 1 atm is
lowered. Surface waves have been implicated as responsible increased by the hydrostatic pressure above the bubble. pA is
for the destruction. As the amplitude of these oscillations can the maximal ultrasound driving amplitude with frequency f
increase to the size of the radial oscillation amplitude, the = 20 kHz. The ambient bubble radius R0 is calculated at pA
bubble splits into fragments 关4兴 and may not survive this = 0. ␴ is the surface tension and ␩ is the viscosity of the
catastrophic event. liquid. The fluid parameters are taken from tabulated values
Radial and surface oscillations of bubbles are calculated 关7兴.
and their type and behavior is determined. The total param- b is a van der Waals hard-core radius and ␥ a polytropic
eter space of sonoluminescing bubbles is scanned for differ- exponent. Its value is set between 1 共=isothermal兲 and the
ent instabilities and the results are mapped. adiabatic exponent of the gas according to an instantaneous
Péclet number 关8–10兴 Pe= R20兩Ṙ共t兲兩 / R共t兲␬, reflecting thermal
II. RADIAL OSCILLATIONS conduction at the involved time scales. ␬ is the thermal dif-
fusivity of the gas. To avoid a change to isothermal behavior
The Gilmore model 关5兴 describing the radial motion of a
at the bubble wall turning point during maximum compres-
bubble in a compressible liquid is integrated numerically.
sion, the bubble wall velocity is kept to its maximum during

冉 冊 冉 冊 冉 冊 冉 冊
RR̈ 1 −
Ṙ
C
3
+ Ṙ2 1 −
2
Ṙ
3C
=H 1+
Ṙ
C
R
+ Ḣ 1 −
C
Ṙ
C
,
positive bubble accelerations at collapse. It is smoothly ap-
proaching the real velocity during the rest of the cycle.
b and ␥ are updated during calculations to reflect the ac-
共1兲 tual gas content. The value of the thermal diffusivity ␬ is
variable, as it depends on the varying density ␳g of the gas:

H= 冕 p共R,R˙兲

p⬁
␳−1
l dp,
p+B
p0 + B
= 再冎
␳l
␳0
n
, 共2兲 ␬=
k
=
␳g共R0兲
␳gc p ␳g共R兲
⫻ 2 ⫻ 10−5 m2 s−1 . 共4兲

k is the thermal conductivity and c p the specific heat at con-

兩C = c兩r=R = 冏冑 冏 冉dp
d␳l r=R
= c0
p共R,Ṙ兲 + B
p0 + B
冊 共n−1兲/2n

,
stant pressure. The value of ␬ is scaled by the ratio of the
ambient gas density to actual density.
The density in the bubble is calculated by

1539-3755/2008/77共6兲/066309共10兲 066309-1 ©2008 The American Physical Society

JOACHIM HOLZFUSS PHYSICAL REVIEW E 77, 066309 共2008兲

ntotalM̄ Ṙ − vH2O
␳g共R兲 = 4 , 共5兲 ⍀= . 共11兲
3 ␲R
3
cpeak

4 3 ntotalRgasTB 4 3 ⍀ is a ratio of velocities, Ṙ is the bubble wall velocity, vH2O

␲R = + ␲b . 共6兲 is the vapor velocity, and cpeak is the velocity belonging to
3 pg 3
the peak of the Maxwell-Boltzmann velocity distribution,

冑
pg is the pressure in the bubble, and Rgas
= 8.3143 J mol−1 K−1 is the gas constant. The van der Waals 2RgasTs
cpeak共Ts兲 = . 共12兲
hard-core parameter b is calculated as an average from the M H2O
tabulated hard-core values and the number of moles of the i
different gases, i = N2, O2, Ar, and H2O, M̄ is the molar mass The change of mass per unit time and unit area j can be
averaged over all gases, and ntotal the total number of moles. expressed as
b and M̄ are updated continuously.
ṅH2OM H2O
The temperature TB is taken to be uniform within the j= = ␳g,H2O共Ṙ − vH2O兲 共13兲
bubble. It is calculated via the adiabatic compression of a 4␲R2
van der Waals gas by
and when inserted into 共11兲 leads to
TB = T0 冉 R30 − b3
R3 − b3
冊 ␥−1
共7兲
⍀=
ṅH2O R
共14兲
nH2O 3cpeak
with the ambient liquid temperature T0. Other more complex
approaches exist that include thermal conduction and a tem- for a spherical bubble volume. For small values of ⍀ the
perature jump across the bubble interface 关11,12兴. approximation ⌫ = 1 − ⍀冑␲ can be made 关23兴, leading to an
Evaporation and condensation of water molecules at the effective constant evaporation coefficient of ␣eff = 2␣ / 共2
bubble wall 关3,12–15兴 are included in the model for the − ␣兲 in 共8兲 and an effective correction factor ⌫eff of unity.
bubble dynamics, as experimental results 关16,17兴 stress the Calculations show, that ⌫ varies by as much as 20% around
importance of a decrease of the polytropic exponent induced unity during the collapses. However, in the observed param-
by water vapor at bubble collapse. A simple Hertz-Knudsen eter range almost no notable difference exists in the amount
model 关18,19兴 for the change of moles nH2O of water vapor in of water vapor at collapse time between results of calcula-
the bubble is tions using a constant value of ␣eff and setting ⌫eff = 1, and of
those using 共8兲 with the variable expressions 共10兲 and 共14兲
4␲R2 ␣c̄共Ts兲 sat together with a smaller value of the evaporation coefficient
evap
ṅH2O = ṅH cond
O − ṅH O = 共␳g,H O − ⌫␳g,H2O兲.
2 2 M H2O 4 2
␣ = 2␣eff / 共2 + ␣eff兲. Therefore a constant value of ␣eff = 0.4
共8兲 关24兴 is taken in the following calculations.

␣ is the constant evaporation coefficient 共also called the ac- III. SURFACE WAVES
commodation coefficient or sticking probability兲,

冑
Not only is a bubble in a liquid oscillating radially, but the
8RgasTs surface may also exhibit a wavy oscillation, as has been ob-
c̄共Ts兲 = 共9兲 served earlier 关25兴. The spherical stability of a bubble has
␲ M H2O
been analyzed by a great number of authors 关26–30兴. The
is the average velocity of molecules of a Maxwell- problem is involved, since during one oscillation the direc-
Boltzmann distribution, ␳g,H2O is the density of water vapor tion of bubble acceleration changes as well as the density
of molar weight M H2O in the bubble, and ␳g,H
sat
is the satu- within the bubble.
2O
rated vapor density 关7兴. The bubble surface temperature is A linear analysis is made studying the perturbation of the
taken as Ts = T0. The density of water vapor ␳g,H2O depends bubble shape 关26,31兴
on the bubble dynamics and is calculated along with the
r共t,⌰, ␾兲 = R共t兲 + an共t兲Y m
n 共⌰, ␾兲, 共15兲
bubble equation with the help of 共5兲. The simple model 共8兲
takes the temperature distributions in the bubble and liquid where R共t兲 is the instantaneous bubble mean radius, Y m n a
as fixed and does not capture all effects occurring during surface harmonic of degree n and order m, and an its ampli-
evaporation, as would more complex treatments 关20–22兴. ⌫ tude. Prosperetti developed a theory for surface oscillation of
is a correction factor for nonequilibrium conditions induced a bubble 关32兴 leading to an integro-differential equation for
by mass motion of vapor and bubble wall movement the amplitude disturbance an. In the linear analysis, the per-
共Schrage correction兲 关11,14,23兴, turbation dynamics is independent of m. For small liquid

⌫ = e−⍀ − ⍀冑␲ 1 −
2
冉 冑冕 2
␲ 0
⍀
2
冊
e−x dx , 共10兲
viscosity integral terms can be dropped. Neglecting the den-
sity and viscosity of the gas content in a liquid with constant
density ␳l = ␳0, a simpler expression is 关33兴

066309-2
SURFACE-WAVE INSTABILITIES, PERIOD DOUBLING,… PHYSICAL REVIEW E 77, 066309 共2008兲

冉
än + 3
Ṙ
R
␩
+ 2共n + 2兲共2n + 1兲 2 ȧn
␳ lR
冊 A. Rayleigh-Taylor Instability
Basically two types of instability have been identified.
The Rayleigh-Taylor 共RT兲 instability 关39,40兴 acts on very

+ 共n − 1兲 − 冉 R̈
R
␴ ␩Ṙ
+ 共n + 1兲共n + 2兲 3 + 2共n + 2兲 3 an
␳ lR ␳ lR
冊 short time scales on 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 accel-
= 0. 共16兲 erated into the heavier fluid, namely, water. Small distur-
bances an Ⰶ R may then be amplified exponentially, eventu-
For a bubble of fixed radius R, the time derivatives vanish
ally leading to a bubble breakoff. This is most prominently
and the natural frequency wn and the damping coefficient bn
taking place around the main bubble collapse. Arguments
for the nth mode can easily be deduced from 共16兲 as
exist in the current literature that the RT instability may be
␴ suppressed by a proper account of thermal effects in the gas
␻2n = 共n − 1兲共n + 1兲共n + 2兲 , 共17兲 关28,33兴. A suppression when comparing to results in 关34兴 is
␳ lR 3
also calculated 关37兴, when the increasing density of the bub-
␩ ble’s interior during the collapse is accounted for 关Eqs.
␤n = 共n + 2兲共2n + 1兲 . 共18兲 共19兲–共22兲兴. Different methods for the calculation of the RT
␳ lR 2 instability have been developed. Sudden growth of the sur-
The first natural frequencies ␻n / 2␲ , n = 2 , 3 , 4, can be cal- face oscillation amplitude to the size of the bubble radius
culated with the usual parameters to 2.638, 4.817, and 7.226 defines the instability threshold, resulting in a breaking of the
MHz for a bubble of 5 ␮m ambient radius and 3.687, 6.731, bubble. To induce this instability numerically, different strat-
and 10.097 MHz for a 4 ␮m bubble. This suggests that the egies exist. In 关34兴 a random displacement of size 0.1 nm is
n = 2 surface mode is the first excited mode due to its lower added as a microscopic fluctuation to an in 共19兲 during inte-
natural frequency. gration. Others have added a fixed initial kick to the surface
Lohse and co-workers 关34,35兴 went one step beyond the oscillation amplitude of 10 nm shortly before the bubble col-
above approximation and arrived at a boundary-layer-type lapse 关29,37兴, added Gaussian distributed noise with 0.1 nm
approximation 共BLA兲, where the vorticity only acts within a standard deviation as molecular fluctuations 关41兴, or refor-
small distance from the bubble. Using their BLA approxima- mulated the problem as a stochastic differential equation
tion on the original equations 关32兴, including the dependence 关36兴. Here, Gaussian distributed noise with 0.5 nm standard
on the time varying gas density ␳g, 关32兴 and neglecting gas deviation is added to 共19兲 during integration. An important
viscosity 关36,37兴, one gets quantity is the absolute value of an共t兲 / R共t兲: A maximum ex-
ceeding a value of unity implies a splitting of the bubble.
än + Bnȧn + Anan = 0, 共19兲 Calculation of the above ratio shows that it depends on the
amount and type of noise. Changing the strength of the fluc-

Bn = 3
Ṙ
R
+
␳l
冉+
n+1 n
␳g
冊 −1 tuation 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
2␩ n共n + 2兲2 1 disturbances induced by the RT instability is inversely pro-
⫻ − 共n − 1兲共n + 2兲 , 共20兲
R 2
n + 1 1 + 2␦/R portional to their wavelength, higher order surface harmonics
are used for the calculation. As a side effect, the afterbounce

An = 冉 ␳l
+
n+1 n
␳g
冊 冋冉
−1
n+2
n
␳g −
n−1
n+1
␳l
R̈
R
冊 instability 共see next paragraph兲 as a parametric effect is dis-
appearing for high-order modes 共Fig. 1兲. The Gilmore equa-
tions together with 共19兲–共22兲 are used in Fig. 1 for an argon-

+ 共n − 1兲共n + 2兲
␴
R
Ṙ
冉
3 + 2␩ 3 共n − 1兲共n + 2兲
R
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

−
n共n − 1兲共n + 2兲
n+1
1
1 + 2␦/R
冊册 , 共21兲
of the amplitude of the n = 6 mode to a value larger than the
threshold of breakup is seen during main bubble collapse.
Parameter space calculations show that the existence of
the RT instability is limited to a small region in phase space
where ␦ is the boundary layer thickness, suggested 关38兴 to be
共Fig. 2兲. The radius for an undriven bubble is taken as the
defined as
equilibrium radius. Due to usage of the Gilmore model,

␦ = min 冉冑 ␩ R
,
␳l␻ 2n
冊 共22兲
which includes a better modeling for acoustic damping when
compared to simpler models, the inclusion of the density
correction in 共19兲–共22兲, and water evaporation-condensation,
where ␻ / 2␲ is the ultrasound driving frequency. In the fol- it is seen that the RT-unstable region is very small and con-
lowing calculations the boundary layer is set to a value ␦ fined to high pressures and small mean equilibrium bubble
= 0, which, as can be seen by the results, models stability radii. The region gets smaller at increasing pressures 共com-
boundaries in fairly good agreement with experimentally pare 关34,36,37,42兴兲. Furthermore, the amount of added noise
measured boundaries 关37兴. in the calculation has to be increased by a factor of 5 to show

066309-3
JOACHIM HOLZFUSS PHYSICAL REVIEW E 77, 066309 共2008兲
Surface Distortion/Mean Radius
1.5

Surface Distortion/Mean Radius

5
n=2
1 0

0.5 -5

0 3
2
1.5 1
n=6 0
1
-1
20 30 40 50
0.5 Time [µs]
0
FIG. 3. Two types of afterbounce instabilities: Upper: Driving
0 0.5 1 1.5 pressure is 1.4 bars, ambient radius 5.23 ␮m. Lower: Driving pres-
Time [µs] sure is 1.7 bars, ambient radius 2.51 ␮m. Driving frequency is 20
kHz. The times of the main collapses are indicated by arrows.
FIG. 1. Rayleigh-Taylor instability: The driving pressure is 1.55
bars, the mean ambient bubble radius 0.903 ␮m, and the driving
frequency 20 kHz. Shown is the ratio of the amplitude of surface may induce very large surface oscillations leading to the de-
distortion to mean radius. Maxima and minima are seen around the struction of a bubble within a single cycle. The instability is
main bubble collapse and subsequent afterbounces. The n = 2 共upper most predominant for low-order surface harmonics. Two
graph兲 surface harmonic shows a smaller distortion at the main types of this instability can be observed 共Fig. 3兲: It can be
collapse which is amplified by the afterbounces. The n = 6 共lower triggered by an initial RT instability 共Fig. 3, lower graph兲
graph兲 surface harmonic shows a large distortion at the main col- which induces a large surface deformation not leading to
lapse and hardly any reaction to the afterbounces. breakup during the main collapse. The deformation is not
damped out until the next afterbounce, thus leading to an
a notable effect. It may be asked if the added amount of 0.5 even larger surface deformation during the first afterbounce
nm average Gaussian noise amplitude is too large for experi- collapse while giving the same short-time behavior as an RT
mental conditions. If the answer is yes then a destruction of instability. This may repeat until the afterbounces finally get
the bubble by a genuine Rayleigh-Taylor instability is not too small for further amplification. This type of instability
happening. can only be demonstrated numerically with a large amount of
noise fluctuations of 0.5 nm standard deviation of Gaussian
distributed noise.
B. Instability by secondary collapses Another type of amplification by afterbounces is seen to
An amplification of the surface mode oscillations by re- start to exist near the bifurcation to parametric instability
peated secondary collapses 共afterbounces兲 subsequent to the 共Fig. 3, upper graph兲: Due to the added noise the surface
main collapse has been termed afterbounce instability 关34兴. It oscillations get very large and may exceed the minimum of
the mean bubble radius thus leading to destruction. This hap-
pens only at later afterbounces which resonantly amplify the
3
surface deformations. These deformations grow only from
10
cycle to cycle once the parameters are set to values above the
Osc. Amplitude Ratio

2.5
1
threshold for parametric instability. This type of afterbounce
instability also exists for small amounts of added noise. Each
Radius [µm]

2 instability exists in its own distinct region in the parameter

0.1

space 共Fig. 4兲.

1.5 0.01

C. Parametric instability
1
The parametric instability, in planar geometry also known
as the appearance of Faraday waves 关43,44兴, describes the
1.3 1.35 1.4 1.45 1.5 1.55 1.6 1.65 1.7 growth of the absolute magnitude of the surface deformation
Pressure [bar]
amplitude an over many cycles of the driving sound. If small
FIG. 2. 共Color online兲 Rayleigh-Taylor instability of an argon- disturbances an Ⰶ R are vanishing from period to period, the
water vapor bubble: In the parameter space spanned by driving bubble is said to be parametrically stable. Some debate has
pressure and equilibrium radius, the average maximum of the ratio occurred concerning the role of viscosity concerning the
of instantaneous surface deformation to mean radius is color coded. parametric instability 关45,46兴 resulting in an underestimation
The maximum deformation ratio of modes n = 2–6 is taken around of the excitation threshold of bubble radii for low driving
the main collapse and averaged over ten consecutive driving oscil- pressures. A semianalytical approach 关29兴 and more exact
lation cycles. Within the contour lines the value is larger than unity, numerical models 关47,48兴 finally showed the validity of Eqs.
implying splitting of the bubble. The equilibrium radius is taken as 共19兲–共22兲 at least in the parameter region for sonolumines-
the value of radius with no periodic forcing amplitude applied. cence.

066309-4
SURFACE-WAVE INSTABILITIES, PERIOD DOUBLING,… PHYSICAL REVIEW E 77, 066309 共2008兲

8 80

10
7

Osc. Amplitude Ratio

60
6

Radius [µm]
1
Radius [µm]

40
5
0.1

4 20
0.01
3
0
2
0 100 200
1 (a) Time [µs]

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

Pressure [bar]

Radius [µm]
60

FIG. 4. 共Color online兲 Amplification of surface mode oscilla- 40

tions of an argon-water vapor bubble by afterbounces: In the pa- 20
rameter space spanned by driving pressure and equilibrium radius

Surf.Amp. / Radius
the average maximum of the ratio of surface deformation 共n = 2兲 0

amplitude to mean radius is color coded. The average maximum 1

occurring during a whole driving oscillation period of 6 different 0
cycles is taken. The individual cycles are nonconsecutive and ran- -1
domly disturbed. Within the contour lines the value is larger than
100 120 140 160 180 200
unity, implying splitting of the bubble. Two main regions of this (b) Time [µs]
instability can be seen to exist.
FIG. 5. Parametric instability of an argon-water vapor bubble.
The parametric instability is calculated in principle using The maximal amplitude of surface deformation increases during
Floquet’s theorem by calculating the amplification rate of a several periods of the driving sound to a value larger than the ra-
small perturbation after one period of a steady-state radial dius, thus leading to bubble breakup. Upper graph: The upper trace
oscillation 关34兴. In the parameter region analyzed here, this shows bubble radius, while the lower trace shows surface oscilla-
can be relaxed to monitoring the growth of the surface de- tion amplitude 共n = 2兲. Equilibrium radius is R0 = 5.74 ␮m, 1.4 bars
formation amplitude an over several cycles. It is found that is the sound pressure. Lower graph: Excerpt of upper graph around
the n = 2 mode has the lowest threshold for a per period am- point of splitoff. f = 20 kHz is the driving frequency for all
plification rate larger than unity. Figure 5 shows a typical calculations.
example of parametric amplification of an n = 2 mode. Dur-
ing a single period the time trace shows the development of
a surface oscillation. It is seen that the surface oscillations within adjacent cells of Fig. 6. Time sections of blowups of
are driven by the afterbounces. The oscillation period equals subsequent oscillations are shown. Both rows show a para-
one-half of the driving afterbounces. The envelope shows an metric amplification of the surface oscillation amplitude.
increase of the oscillation amplitude followed by a subse- However, in contrast to the second, the first row shows an
quent decrease. The remaining amplitude of the surface de- alternating sign symmetry of the wave form of consecutive
formation at the time of the next main collapse is amplified periods. This also results in an opposite sign of the amplitude
by the RT instability, again resulting in a parametric oscilla-
tion. If the maximal amplitude is increasing over the time 8
scale of several periods this leads to bubble breakup. 100
Figure 6 shows the value of the growth rate per driving 7 1
Growth Rate

0.01
period of the n = 2 surface oscillation mode in the parameter 1x10 -4
Radius [µm]

space of SBSL. Parameters of the equilibrium radius and 6 1x10 -6

1x10 -8
driving pressure from within or above the white contour lines 1x10 -10
show parametric instability leading to bubble breakup. As a 5 1x10 -12

tendency, a lower driving amplitude and smaller bubbles are

more stable. However, a fine structure is superimposed onto 4
the graph, showing individual cells separated by horizontal
and vertical lines. At the lines the amplification rates are 3
smaller compared to those within each cell. This is true be- 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
Pressure [bar]
low the stability threshold as well as above it. Near the
boundary a more complicated scenario exists: Eventually a FIG. 6. 共Color online兲 Parametric stability of a argon-water va-
cell contains an island surrounded by a white contour line por bubble: In the parameter space spanned by driving pressure and
denoting the presence of an unstable region within an other- equilibrium radius the growth rate of surface deformation per pe-
wise stable cell. riod 共n = 2兲 is color coded. Within the white contour lines the value
Figure 7 shows simulations of the amplitude of the n = 2 is larger than unity, implying parametric growth of the surface
surface oscillation mode with two parameter settings from deformation.

066309-5
JOACHIM HOLZFUSS PHYSICAL REVIEW E 77, 066309 共2008兲

also its spatial position and driving pressure change, a bubble

Surface Deformation / Radius [x10 ]
3
-3

2 may show different radiation behavior, changing from uni-

1 form to nonuniform, and vice versa. This may be understood
0 as radiation from a bubble entering and leaving islands of
-1 symmetrically different parametric stability. Interestingly,
-2 measurements in 关51,52兴 show a double periodic cycle in the
-3
0 10 20 30 40 60 70 80 90 110 120 130 140 160 170 180 190
intensity of spatially directed radiation near the upper thresh-
3 old of SBSL and during unstable SBSL. This can be well
2 understood with the results of Fig. 7, where the surface am-
1 plitude at collapse alternates for subsequent periods of the
0 driving, revealing a symmetry change of the bubble surface
-1 at main collapse with a doubled period. As the average radial
-2 oscillation of a bubble is only driving the surface oscillation
-3 and experiences no back coupling, it would be unaffected in
0 10 20 30 40 60 70 80 90 110 120 130 140 160 170 180 190
this scenario, resulting in a constant periodic collapse as
Time [µs] measured in 关52兴. The absolute value of the parametric sta-
bility is only slightly above one in the first two 共for increas-
FIG. 7. Parametric growth of surface deformation of an argon-
water vapor bubble: The time dependence of the 共n = 2兲 surface
ing radius兲 elongated cells above the threshold in Fig. 6.
oscillation per radius ratio is shown for four consecutive time sec- Thus it may be possible that nonlinear surface mode effects
tions around the main collapse 共marked by arrows兲. Data from two may inhibit the growth, leading to bubble destruction. Also
different parameter settings have been taken. First row: equilibrium in the case of a bubble growing by diffusion 关4兴 it may pass
radius R0 = 5.839 ␮m at 1.35 bars; second row: R0 = 6.0 ␮m at 1.35 through an unstable island fast enough without being de-
bars. Note the alternating symmetry in the first row. stroyed 关52,53兴.
The existence of the cells in Fig. 6 is a consequence of
coupling of nonlinear oscillators 共one sinusoidally driven ra-
during the main collapse 共marked by arrows兲. Whenever a
dial oscillator driving the surface oscillator兲 and does not
parameter setting is picked from an adjacent cell, a transition
depend on the types of models and approximations used.
from nonalternating to alternating behavior is made. This is
Signs of a cell-like structure can also be seen in graphs
true within the whole parameter space, for parametrically
showing parametric boundaries published by other authors
stable or unstable settings. An alternating parametrically am-
using different models and approximations 关27,33,38兴. The
plified wave form is also seen in Fig. 5.
exact positions of lines of instability in parameter space
The amplitude of the n = 2 oscillation is about 2% of the
spanned by ambient radius and driving pressure do depend
mean radius at main collapse at the time the ratio of n = 2
somewhat on the models and parameters used 共in accordance
oscillation to mean radius exceeds 1 during the afterbounce
with the current literature兲. Their selection can only be jus-
oscillations, implying breakup of the bubble.
tified by comparison with experimental data. Models using a
nonuniform bubble interior should find inner high-pressure
IV. DISCUSSION peaks or shock waves being detached from the elliptical
The region of hydrodynamic stability of driven bubbles in bubble surface, resulting in a nonspherical core. This will
the parameter range of sonoluminescence has boundaries in lead to a nonisotropic emission as well. In fact, molecular
the direction of increasing equilibrium radius and increasing dynamics simulations have been published 关54兴 that report
driving pressure. The Rayleigh-Taylor instability is numeri- nonisotropic emission upon collapse of slightly ellipsoid
cally shown to be limited to a region of very small bubbles at bubbles during SBSL.
high driving pressures. The instability set up by secondary An instability that has as yet not been analyzed in SBSL
collapses limits the accessible parameter space also in the research, but may be of importance, is the Richtmyer-
direction of larger equilibrium radii as well as higher pres- Meshkov instability 关55–57兴. It occurs when a shock wave
sures. The determination of the exact position of the bound- interacts with the gas-fluid interface. The Kelvin-Helmholtz
aries is hindered by the amount of added noise used in the instability 关57兴 occurs when there is a fluid motion parallel to
numerical procedure. When comparing to experimentally de- the interface between gas and water, which may be present in
termined boundaries, the experimentally existing noise the form of microstreaming 关58,59兴.
should also be considered.
The parametric instability limits the parameter space V. APPROXIMATE UNIVERSAL BOUNDARY
mostly in the direction of growing radii. A consequence for a
sonoluminescing bubble driven near the threshold of para- The calculated threshold for the parametric instability us-
metric stability is that an appreciable distortion 共e.g., prolate- ing the boundary layer approximation in 共22兲 is found to be
oblate ellipsoid兲 exists during the collapse. This may result in lower than the threshold at which a bubble disappears in
a directional radiation of sonoluminescence light output 关49兴. experiments. It has been shown 关37兴 that smaller values of
This has been seen previously in experiments 关50–52兴. Mea- the size of the boundary layer ␦ 关Eq. 共22兲兴 result in a shift of
surements presented in 关51–53兴 have shown that during un- the parametric instability to larger ambient radii. Actually,
stable SBSL, where the equilibrium radius of the bubble and setting ␦ = 0 resulted in the best coincidence of numerical and

066309-6
SURFACE-WAVE INSTABILITIES, PERIOD DOUBLING,… PHYSICAL REVIEW E 77, 066309 共2008兲

2.5
RPImax f = -2.296+2.648*(pA/p0)
2 1.5
RPImax f [m/s]

RPImax f [m/s]
1.5
1
1
Krefting 23 kHz, 1 atm + 230 Pa
Miniforscher 24 kHz
0.5 Thomas 14.62 kHz
0.5 Vazquez 33.8 kHz, 2atm
Barber 26.4 kHz, 1 atm
0 Gaitan 20.6 kHz, 1 atm + 600 Pa
1 1.2 1.4 1.6 1.8 Levinsen 21.93 kHz
pA / p0
0
1.1 1.2 1.3 1.4 1.5
FIG. 8. Maximal bubble radius times driving frequency at para- pA / p0
metric instability threshold as a function of ratio of driving pressure
over ambient pressure. Bubble oscillations below the line are para- FIG. 9. Maximal bubble radius times driving frequency as a
metrically stable. Calculations are shown for ambient pressures of function of normalized pressure for various data sets obtained from
+ / −1.4 kPa Ar 共squares, triangles兲 around 1 atm 共circles兲 with a literature: crosses 关60,61兴, stars 关62兴, filled circles 关63兴, down tri-
driving frequency of 20 共filled symbols兲, 14 共hollow squares兲, 24 angles 关16兴, up triangles 关64兴, filled squares, own data, hollow
共stars兲, and 30 kHz 共hollow circles兲 at 1 atm. The boundary layer is circles 关65兴. Where indicated the ambient pressure 共set to 1 atm if
␦ = 0. The line shows a linear regression over all data points. unknown兲 has been augmented by hydrostatic pressure guessed
from published cell dimensions. The line is determined from Fig. 8.
experimental stability boundaries. One has to cope with the
fact that it is somewhat difficult to measure the bubble’s
pressure 关67,68兴. The authors find a linear increase of the
ambient radius. This is because of limited optical resolution
number of photons as the ambient pressure p0 is decreased.
and the question of when and how to measure an ambient
radius in the case of sinusoidal driving. A solution is to mea- Here we are concerned about the maximally obtainable
sure the maximal radius Rmax of the bubble during an oscil- bubble radius close to the border of instability. At constant
lation. frequency and driving pressure, a decrease of p0 increases
PI
Figure 8 shows Rmax f, which is the maximal bubble radius the maximal bubble radius at the border of parametric insta-
at the first lower border of parametric instability times the bility too, eventually leading to a higher number of emitted
driving frequency, as a function of the driving pressure over photons.
ambient pressure ratio. An almost linear dependence is sug- The nearly perfect linear dependence of the threshold line
gested over a large parameter variation covering the whole
PI
Rmax f as a function of normalized pressure pA / p0 in the pa-
SBSL driving parameter space. A linear regression over all rameter range for SBSL suggested by Figs. 8 and 9 can be
data points results in RmaxPI
f = −2.296+ 2.648共pA / p0兲 with a explained: A larger maximal radius leads to a more violent
standard error of 0.01 in both constants. A more detailed bubble collapse. The resulting larger amplitudes of after-
calculation 共Fig. 11兲 of the boundary shows the same but bounces lead to stronger excitation of surface waves. The
very squeezed roughness as Fig. 6, with all the long island increasing frequency of the collapses leads to a further am-
lines running parallel to the line in Fig. 8. plification as the surface waves have less time to be damped
Various data of the maximal bubble radius at the border of between each excitation. The dependent axis has units of
SBSL extinction have been obtained from the literature and velocity. Parametric instability is induced when this thresh-
are shown in Fig. 9. The data sets have been taken by differ- old velocity is passed over. The requirement of pressure axis
ent authors, with different setups, driving frequencies, and normalization by the ambient pressure is also seen in the
pressures. Ambient pressures have not always been pub- case of the ambient pressure dependence of diffusive stabil-
lished, nor has the hydrostatic pressure induced by the water ity lines 关3,69,70兴. A larger ambient pressure decreases the
column on top of the bubble, which has hence been set to an maximum radius during the oscillation.
elaborated estimate. A very good agreement with the calcu- The existence of the approximate universal boundary of
lated linear relationship in Fig. 9 is seen. In particular, the the parametric instability is shown to be robust against
detailed results from 关63兴 agree well. The outmost lying changes of the bubble equation and various approximations
circle corresponds to a published ambient pressure being made. Figure 10 shows parametric instability lines for the
somewhat “out of sequence.” The down triangles 共from re- Gilmore, Keller-Miksis 关71兴, and modified Rayleigh-Plesset
sults of Fig. 5, curves a and b in 关16兴兲 lie outside because equations 关3兴 driven at f = 20 kHz. The equations are aug-
they are taken at an ambient pressure of 2 bars, which is not mented with the other model equations mentioned earlier.
in the applicability range of the dependence suggested by Hardly any difference with the choice of bubble equation is
Fig. 8. The results represented by the squares stem from a visible. Due to its more complete modeling capabilities for
rather noisy experiment. The hollow circle is a maximally radiated sound, the Gilmore equation shows slightly smaller
driven nitrogen bubble 关65,66兴 at 9 °C close to breakup. values of the maximal radius for large forcing amplitudes
Experimental and numerical studies have been made con- due to acoustic energy radiation upon collapse in the form of
cerning the dependence of light output on change of ambient shock waves. Calculations for the modified Rayleigh-Plesset

066309-7
JOACHIM HOLZFUSS PHYSICAL REVIEW E 77, 066309 共2008兲

2.5 100

PM AB1

Maximal Radius [µm]

2 80 AB2
RPImax f [m/s]

1.5 60

RPm adiabatic
1 RPm BLA 40
RPm ideal gas
RPm noevap
RT
KM
0.5 RPm 20
Gilmore
fit
0 0
1.1 1.2 1.3 1.4 1.5 1.6 1.7 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
pA / p0 Pressure [bar]

FIG. 11. 共Color online兲 Phase space spanned by driving pressure

FIG. 10. Maximal bubble radius times driving frequency at
and maximal radius showing regions of Rayleigh-Taylor 共RT兲, af-
parametric instability threshold as a function of normalized pres-
terbounce 共AB兲, and parametric 共PM兲 instabilities. A region for high
sure. Data from various models and approximations are shown.
forcing amplitudes exists between AB1 and AB2 where none of the
Circles are calculated from the Gilmore bubble model using all
above instabilities exist.
equations specified earlier. Crosses denote the Keller-Miksis model
共KM兲. The modified Rayleigh-Plesset bubble model 共RPm兲 is Levinsen et al. In order to find an explanation for this, dif-
shown with several changes applied to demonstrate the robustness ferent hydrodynamical instabilities have been mapped
of the approximation of a universal boundary. The line is deter- throughout the whole parameter space. The parametric insta-
mined from Fig. 8. bility induces slowly growing surface waves, calculated as a
prolate-oblate oscillation. The structure of the growth rate of
model with an additional change in approximation have been this instability mapped in parameter space shows a ragged
made: A purely adiabatic bubble compression with constant threshold line, consisting of multiple islands separated by
adiabatic exponent 关Eq. 共7兲兴 shows the largest deviation from stability zones. The dynamical behavior in adjacent islands is
the line calculated from Fig. 8 to values that are also well different. In one island the surface oscillations catch up with
below the experimental data points from Fig. 9. The neglect each other such that the bubble is always prolate or always
of the evaporation-condensation equations also results in oblate at collapse. In the neighboring island, a period-2 os-
overall smaller values, as does the inclusion of the BLA cillation is seen such that for one period the bubble is oblate
approximation via Eq. 共22兲. The latter is in accordance with while in the next it is prolate at collapse. This last type of
oscillation will result in a spatially directed emission of light
the results from 关37兴 on the ambient radius. Results for a
with a period-doubled intensity, just as observed in the ex-
different equation of state for the gas inside the bubble 共ideal
periment.
gas兲 do not show any considerable difference from the re- The Rayleigh-Taylor instability occurring around collapse
sults for a van der Waals hard-core gas. This motivates the upon acceleration of the lighter fluid 共gas兲 of the bubble into
assumption that using a more realistic, but computationally the heavier fluid 共water兲 is calculated using varying gas den-
expensive, soft shell model 关72兴 also leads to a very small sity and higher-order disturbances. This instability is found
difference between the results. to occur only at small ambient radii and stronger forcing, and
Simulations using molecular dynamics have been pub- should therefore not occur at the upper driving pressure
lished 关54,73,74兴 that use a bubble equation as a launch con- threshold of SBSL. Two different types of afterbounce insta-
dition. Back reactions of nonisotropic bubble interiors on bility have been identified and shown to exist in different
surface oscillations have not yet been studied in this context. regions of the parameter space. For experiments with a high
Analysis of those may lead to further agreement of theory amount of noise 共sound disturbances, higher ambient tem-
and experiment. perature兲 one type of afterbounce instability may impose an
Figure 11 shows the regions of all instabilities analyzed upper limit to SBSL.
before in a single graph. The comparison shows that the A global threshold function has been found in the driving
threshold for the noise-induced afterbounce instability AB1 pressure–maximum radius space at the parametric instability.
is shifted to smaller maximal radii below the parametric in- Within the limits of the ragged islands mentioned earlier,
stability. The Rayleigh-Taylor instability RT occurs at which are visible only in very detailed calculations, a rescal-
smaller maximum radii near the border of the afterbounce ing of the pressure axis with the ambient pressure and the
instability AB2. A small gap with none of the above insta- maximum radius at threshold with the driving frequency re-
bilities extends to larger maximal radii with increased driv- veals a straight line dependence in the SBSL parameter
ing pressure. range. All available experimental data from the literature
have been inserted into a graph. Within experimental accu-
racy a global threshold value of the parametric instability as
a function of rescaled pressure can be deduced.
VI. CONCLUSIONS
ACKNOWLEDGMENT
Nonisotropic radiation in combination with spatially di-
rected period doubling has been observed in experiments by The author thanks M. Levinsen.

066309-8
SURFACE-WAVE INSTABILITIES, PERIOD DOUBLING,… PHYSICAL REVIEW E 77, 066309 共2008兲

关1兴 D. F. Gaitan, L. A. Crum, C. C. Church, and R. A. Roy, J. 2808 共1996兲.

Acoust. Soc. Am. 91, 3166 共1992兲. 关35兴 S. Hilgenfeldt, D. Lohse, and M. P. Brenner, Phys. Fluids 9,
关2兴 B. P. Barber, R. A. Hiller, R. Löfstedt, S. J. Putterman, and K. 2462 共1997兲.
R. Weninger, Phys. Rep. 281, 65 共1997兲. 关36兴 U. H. Augsdörfer, A. K. Evans, and D. P. Oxley, Phys. Rev. E
关3兴 M. P. Brenner, S. Hilgenfeldt, and D. Lohse, Rev. Mod. Phys. 61, 5278 共2000兲.
74, 425 共2002兲. 关37兴 L. Yuan, C. Y. Ho, M.-C. Chu, and P. T. Leung, Phys. Rev. E
关4兴 J. Holzfuss, Phys. Rev. E 71, 026304 共2005兲. 64, 016317 共2001兲.
关5兴 F. R. Gilmore, Hydrodynamics Laboratory, California Institute 关38兴 M. P. Brenner, D. Lohse, and T. F. Dupont, Phys. Rev. Lett.
of Technology, Pasadena, CA, Report No. 26-4, 1952 共unpub- 75, 954 共1995兲.
lished兲. 关39兴 G. I. Taylor, Proc. R. Soc. London, Ser. A 201, 192 共1950兲.
关6兴 J. Holzfuss, M. Rüggeberg, and A. Billo, Phys. Rev. Lett. 81, 关40兴 D. J. Lewis, Proc. R. Soc. London, Ser. A 202, 81 共1950兲.
5434 共1998兲. 关41兴 H. Lin, B. D. Storey, and A. J. Szeri, Phys. Fluids 14, 2925
关7兴 At 20 °C the following parameters have been used in the 共2002兲.
calculations: surface tension ␴ = 0.0725 N m−1, viscosity ␩ 关42兴 H. Lin, B. D. Storey, and A. J. Szeri, J. Fluid Mech. 452, 145
= 0.001 N s m−2, saturated vapor density ␳g,H
sat
共2002兲.
2O
−3
= 0.0173 kg m , diffusion constant DAr = 2.0⫻ 10 m s , −9 2 −1 关43兴 M. Faraday, Philos. Trans. R. Soc. London 121, 299 共1831兲.
diffusion constant DO2 = 1.9⫻ 10−9 m2 s−1, diffusion constant 关44兴 L. Rayleigh, Philos. Mag. 16, 50 共1883兲.
DN2 = 1.9⫻ 10−9 m2 s−1, solubility Ar= 0.061 kg m−3, solubil- 关45兴 S. J. Putterman and P. H. Roberts, Phys. Rev. Lett. 80, 3666
ity O2 = 0.0444 kg m−3, solubility N2 = 0.020 kg m−3, water 共1998兲.
density ␳0 = 998.23 kg m−3. 关46兴 M. P. Brenner, T. F. Dupont, S. Hilgenfeldt, and D. Lohse,
关8兴 D. Lohse and S. Hilgenfeldt, J. Chem. Phys. 107, 6986 共1997兲. Phys. Rev. Lett. 80, 3668 共1998兲.
关9兴 A. Prosperetti, J. Acoust. Soc. Am. 61, 17 共1977兲. 关47兴 P. H. Roberts and C. C. Wu, Phys. Fluids 10, 3227 共1998兲.
关10兴 M. S. Plesset and A. Prosperetti, Annu. Rev. Fluid Mech. 9, 关48兴 C. C. Wu and P. H. Roberts, Phys. Lett. A 250, 131 共1998兲.
145 共1977兲. 关49兴 A. Madrazo, N. García, and M. Nieto-Vesperinas, Phys. Rev.
关11兴 K. Yasui, J. Acoust. Soc. Am. 98, 2772 共1995兲. Lett. 80, 4590 共1998兲.
关12兴 K. Yasui, J. Phys. Soc. Jpn. 66, 2911 共1997兲. 关50兴 K. Weninger, S. J. Putterman, and B. P. Barber, Phys. Rev. E
关13兴 V. Kamath, A. Prosperetti, and F. N. Egolfopoulos, J. Acoust. 54, R2205 共1996兲.
Soc. Am. 94, 248 共1993兲. 关51兴 M. T. Levinsen, N. Weppenaar, J. S. Dam, G. Simon, and M.
关14兴 B. D. Storey and A. J. Szeri, Proc. R. Soc. London, Ser. A 456, Skogstad, Phys. Rev. E 68, 035303共R兲 共2003兲.
1685 共2000兲. 关52兴 J. S. Dam and M. T. Levinsen, Phys. Rev. Lett. 94, 174301
关15兴 B. D. Storey and A. J. Szeri, Phys. Rev. Lett. 88, 074301 共2005兲.
共2002兲. 关53兴 J. S. Dam, M. T. Levinsen, and M. Skogstad, Phys. Rev. E 67,
关16兴 G. E. Vazquez and S. J. Putterman, Phys. Rev. Lett. 85, 3037 026303 共2003兲.
共2000兲. 关54兴 B. Metten and W. Lauterborn, in Nonlinear Acoustics at the
关17兴 R. Toegel, B. Gompf, R. Pecha, and D. Lohse, Phys. Rev. Lett. Turn of the Millennium, ISNA 15, edited by W. Lauterborn and
85, 3165 共2000兲. T. Kurz 共AIP, New York, 2000兲, pp. 429–432.
关18兴 H. Hertz, Ann. Phys. 17, 178 共1882兲. 关55兴 R. D. Richtmyer, Commun. Pure Appl. Math. 13, 297 共1960兲.
关19兴 M. Knudsen, Kinetic Theory of Gases 共Methuen and Co., Lon- 关56兴 E. E. Meshkov, Fluid Dyn. 4, 101 共1969兲.
don, 1950兲. 关57兴 S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability
关20兴 M. Bond and H. Struchtrup, Phys. Rev. E 70, 061605 共2004兲. 共Dover Publications, New York, 1981兲.
关21兴 D. J. E. Harvie and D. F. Fletcher, J. Heat Transfer 123, 486 关58兴 T. G. Leighton, The Acoustic Bubble 共Academic Press, Lon-
共2001兲. don, 1994兲.
关22兴 G. Hauke, D. Fuster, and C. Dopazo, Phys. Rev. E 75, 066310 关59兴 T. Verraes, F. Lepoint-Mullie, T. Lepoint, and M. S. Longuet-
共2007兲. Higgins, J. Acoust. Soc. Am. 108, 117 共2000兲.
关23兴 R. W. Schrage, A Theoretical Study of Interphase Mass Trans- 关60兴 D. Krefting, R. Mettin, and W. Lauterborn, Phys. Rev. Lett.
fer 共Columbia University Press, New York, 1953兲. 91, 174301 共2003兲.
关24兴 I. W. Eames, N. J. Marr, and H. Sabir, Int. J. Heat Mass Trans- 关61兴 The upper limit of the driving pressure is given by 1.34 bars
fer 40, 2963 共1997兲. 关R. Mettin 共private communication兲兴.
关25兴 M. Kornfeld and L. Suvorov, J. Appl. Phys. 15, 495 共1944兲. 关62兴 D. F. Gaitan and R. G. Holt, Phys. Rev. E 59, 5495 共1999兲.
关26兴 H. W. Strube, Acustica 25, 289 共1971兲. 关63兴 C. A. Thomas, R. A. Roy, and R. G. Holt, Phys. Rev. E 70,
关27兴 Y. An, T. Lu, and B. Yang, Phys. Rev. E 71, 026310 共2005兲. 066301 共2004兲.
关28兴 A. Prosperetti and Y. Hao, Philos. Trans. R. Soc. London, Ser. 关64兴 B. P. Barber, C. C. Wu, R. Löfstedt, P. H. Roberts, and S. J.
A 357, 203 共1999兲. Putterman, Phys. Rev. Lett. 72, 1380 共1994兲.
关29兴 V. A. Bogoyavlenskiy, Phys. Rev. E 62, 2158 共2000兲. 关65兴 J. Holzfuss and M. T. Levinsen, Phys. Rev. E 77, 046304
关30兴 A. A. Aganin and T. S. Guseva, J. Appl. Mech. Tech. Phys. 46, 共2008兲.
471 共2005兲. 关66兴 M. T. Levinsen and J. S. Dam, Europhys. Lett. 80, 27004
关31兴 M. S. Plesset, J. Appl. Phys. 25, 96 共1954兲. 共2007兲.
关32兴 A. Prosperetti, Q. Appl. Math. 34, 339 共1977兲. 关67兴 L. Kondic, C. Yuan, and C. K. Chan, Phys. Rev. E 57, R32
关33兴 Y. Hao and A. Prosperetti, Phys. Fluids 11, 1309 共1999兲. 共1998兲.
关34兴 S. Hilgenfeldt, D. Lohse, and M. P. Brenner, Phys. Fluids 8, 关68兴 M. Dan, J. D. N. Cheeke, and L. Kondic, Phys. Rev. Lett. 83,

066309-9
JOACHIM HOLZFUSS PHYSICAL REVIEW E 77, 066309 共2008兲

1870 共1999兲. 关72兴 G. A. Bird, Molecular Gas Dynamics and the Direct Simula-
关69兴 A. Eller and H. G. Flynn, J. Acoust. Soc. Am. 37, 493 共1965兲. tion of Gas Flows 共Oxford University Press, New York, 1998兲.
关70兴 R. Löfstedt, K. R. Weninger, S. Putterman, and B. P. Barber, 关73兴 S. J. Ruuth, S. Putterman, and B. Merriman, Phys. Rev. E 66,
Phys. Rev. E 51, 4400 共1995兲. 036310 共2002兲.
关71兴 J. B. Keller and M. Miksis, J. Acoust. Soc. Am. 68, 628 关74兴 A. Bass, S. Putterman, B. Merriman, and S. J. Ruuth, J. Com-
共1980兲. put. Phys. 227, 2118 共2008兲.

066309-10