Interface dynamics of microscopic cavities in water
Joachim Dzubiella1, ∗
Physics Department, Technical University Munich, 85748 Garching, Germany
(Dated: November 4, 2018)
An analytical description of the interface motion of a collapsing nanometer-sized spherical cavity
in water is presented by a modification of the Rayleigh-Plesset equation in conjunction with ex-
plicit solvent molecular dynamics simulations. Quantitative agreement is found between the two
approaches for the time-dependent cavity radius R(t) at different solvent conditions while in the con-
tinuum picture the solvent viscosity has to be corrected for curvature effects. The typical magnitude
of the interface or collapse velocity is found to be given by the ratio of surface tension and fluid vis-
cosity, v ≃ γ/η, while the curvature correction accelerates collapse dynamics on length scales below
the equilibrium crossover scales (∼1nm). The study offers a starting point for an efficient implicit
modeling of water dynamics in aqueous nanoassembly and protein systems in nonequilibrium.
I. INTRODUCTION
Hydrophobic hydration in equilibrium is a phe-
nomenon which exhibits qualitatively different behavior
at small and large length scales.1,2 While small solutes
(radii R .1nm) are accommodated by water with only
minor perturbations, larger solutes (R &1nm) induce ma-
jor rearrangements of water interfacial structure. As a
consequence the solvation free energy G(R) of small hy-
drophobic cavities scales with solute volume while for
larger cavities it grows with surface area (as a good
approximation near liquid-vapor coexistence) accompa-
nied by weak solvent dewetting at extended restrain-
ing hydrophobic surfaces.3 Furthermore, water, which
is close to the liquid-vapor transition at normal condi-
tions, can minimize interface area by locally evaporat-
ing and forming a ’nanobubble’ within hydrophobic con-
finement. Evidence of bubble formation in confined ge-
ometry has been given early by computer simulations
of smooth plate-like solutes,4 but more recently it has
been demonstrated in varying degrees in atomistically
resolved plate-like solutes,5,6 hydrophobic tubes and ion
channels,7,8 and in the collapse of proteins,9,10 suggesting
that it plays a key role in the stabilization and folding
dynamics of certain classes of biomolecules.11,12 Experi-
mental evidence of nanobubbles in strong confinement (in
contrast to bubbles at a single planar surface3) has been
given for instance in studies of water between hydropho-
bic surfaces,13 in zeolites and silica nanotubes,14,15 and
on a subnanometer scale in nonpolar protein cavities.16
The dewetting induced change in solvation energy is
typically estimated using simple macroscopic arguments
as known from capillarity theory, e.g. by describing in-
terfaces with Laplace-Young (LY) type of equations.14,17
Recently an extension of the LY equation has become
available which extrapolates to microscopic scales by in-
cluding a curvature correction to the interface tension
and considering atomistic dispersion and electrostatic
potentials of the solvated solute explicitly.18 Although
those macroscopic considerations (e.g,. the concept of
surface tension) are supposed to break down on atom-
istic scales they show surprisingly good results for the
solvation energy of microscopic solutes, e.g. alkanes and
noble gases, and quantitatively account for dewetting
effects in nanometer-sized hydrophobic confinement.19
While we conclude that the equilibrium location of the
solute-solvent interface seems to be well described by
those techniques, nothing is known about the interface
dynamics of evolution and relaxation. In this study we
address two fundamental questions: First, what are the
equations which govern the interface motion on atomistic
(∼1nm) scales? Secondly, does the dynamics exhibit any
signatures of the length scale crossover found in equilib-
rium?
On macroscopic scales the collapse dynamics of a (va-
por or gas) bubble is related to the well-known phe-
nomenon of sonoluminescence.20 The governing equa-
tions can be derived from Navier-Stokes and capillarity
theory and are expressed by the Rayleigh-Plesset (RP)
equation.21 We will show that the RP equation simpli-
fies in the limit of microscopic cavities and can be ex-
tended to give a quantitative description of cavity inter-
face dynamics on nanometer length scales. We find a
qualitatively different dynamics than the typical “mean-
curvature flow” description of moving interfaces,22 in par-
ticular a typical magnitude of interface or collapse veloc-
ity given by the ratio of surface tension and fluid viscosity,
v ≃ γ/η. Our study is restricted to the generic case of
the collapse of a spherical cavity and is complemented
by explicit solvent molecular dynamics (MD) computer
simulations. We note here that recently, Lugli and Zer-
betto studied nanobubble collapse in ionic solutions by
MD simulations on similar length scales.23 While their
MD data compares favorably with our results their in-
terpretation and conclusions in terms of the RP equa-
tion are different. We will resume this discussion in the
conclusion section.
In this study we show that a simple analytical approach
quantitatively describes microscopic cavity collapse for
a variety of different solvent situations while the sim-
ulations suggest that the solvent viscosity needs to be
corrected for curvature effects. Our study might offer a
simple starting point for an efficient implicit modeling
of water dynamics in aqueous nanoassembly and protein
http://arxiv.org/abs/0704.0239v1
systems in nonequilibrium.
II. THEORY
The Rayleigh-Plesset equation for the time evolution
of a macroscopic vapor bubble with radius R(t) can be
written as21
= ∆P + 4η
, (1)
where ρm is the solvent mass density, ∆P = P − Pv the
difference in liquid and vapor pressures, η the dynamic
viscosity, and γ the liquid-vapor interface tension. While
for macroscopic bubble radii the inertial terms (left hand
side) control the dynamics, for decreasing radii the fric-
tional and pressure terms (right hand side) grow in rel-
ative magnitude and eventually dominate, so that com-
pletely overdamped dynamics can be assumed on atom-
istic scales:
Ṙ ≃ −
. (2)
A rough estimate for the threshold radius Rt below which
friction dominates is given when the Reynolds number
R = vRρm/η becomes unity and viscous and inertial
forces are balanced. With a typical initial interface ve-
locity of the order of v ∼ γ/η [from R̈(0) = 0 in eq. (1)]
we obtain
Rt = η
2/(ρmγ) (3)
which is ≃ 10nm for water at normal conditions. Note
that this threshold value can deviate considerably for a
fluid different than water and that the viscosity typically
has a strong temperature (T ) dependence which implies
that Rt can change significantly with T .
In equilibrium (Ṙ = 0) the remaining expression in
eq. (2) is the (spherical) LY equation ∆P + 2γ/R = 0.
Thus eq. (2) describes a linear relationship between cap-
illary pressure and interface velocity where R/(4η) plays
the role of an interface mobility (inverse friction).22 In-
terestingly, the mobility is linear in bubble radius which
leads to a constant velocity driven by surface tension in-
dependent of radius (assuming P ≃ 0); this has to be con-
trasted to the typically used capillary dynamics which is
proportional to the local mean curvature ∝ 1/R.22
Generalizations of the LY equation to small scales are
available by adding a Gaussian curvature term (∼ 1/R2)
as shown by Boruvka and Neumann24; that has been
demonstrated to be equivalent to a first order curvature
correction in surface tension, i.e. γ(R) = γ∞(1−δT/R),
where δT is the Tolman length
25 and γ∞ the liquid-vapor
surface tension for a planar interface (R = ∞). The
Tolman length has a magnitude which is usually of the
order of the size of a solvent molecule. Furthermore, it
has been observed experimentally that the viscosity of
strongly confined water can depend on the particular na-
ture of the surface/interface.26 We conclude that in gen-
eral one has to anticipate that - analogous to the surface
tension - the effective interface viscosity obeys a curva-
ture correction in the limit of small cavities due to water
restructuring in the first solvent layers at the hydropho-
bic interface. In the following we make the simple first
order assumption that the correction enters eq. (2) also
linear in curvature (∼ 1/R) yielding
Ṙ = −
∆PR+∆Pδvis + 2γ∞ +
, (4)
where the constant δvis is the coefficient for the first order
curvature correction in viscosity and η∞ the macroscopic
bulk viscosity. Additionally, we define δ = δvis − δT and
second order terms in curvature are neglected. We note
that the choice of the 1/R-scaling of the viscosity curva-
ture correction has no direct physical justification and is
arbitrary. We think however, that a curvature correction
based on an expansion in orders of mean curvature is the
simplest and most natural way for such a choice.
In water at normal conditions the pressure terms in
(4) are negligible so that for large radii (R ≫ δ) the in-
terface velocity is constant and R(t) = R0 − γ∞/(2η∞)t.
This leads to a collapse velocity of about v ≃0.4Å/ps
(40m/s) which is 6% of the thermal velocity of water
vth =
3kBT/m showing that dissipative heating of the
system is relatively weak on these scales. A rough esti-
mate for the dissipation rate can be made by the released
interfacial energy dG(R, t)/dt ≃ d(4πR(t)2γ∞)/dt =
−4πγ2
R(t)/η∞ yielding for instance dG(R, t = 0)/dt ≃
−35kBT/ps for a bubble with R0 =2nm. At small
radii (R ≃ δ) the solution of (4) goes as R(t) ∼
const− (δγ∞/η∞)t decreasing or increasing the ve-
locity depending on the sign of δ = δvis − δT, i.e. the
acceleration depends on the particular sign and mag-
nitude of the curvature corrections to surface tension
and viscosity. For large pressures and radii the first
term dominates which gives rise to an exponential de-
cay R(t) ∼ exp[−∆P/(4η∞)t]. While extending to small
scales we have assumed that the time scale of internal in-
terface dynamics, i.e. hydrogen bond rearrangements,27
is much faster than the one of bubble collapse.
III. MD SIMULATION
In order to quantify our analytical predictions we
complement the theory by MD simulations using ex-
plicit SPC/E water.28 The liquid-vapor surface tension
of SPC/E water has been measured and agrees with the
experimental value for a wide range of temperatures.29
For T = 300K and P = 1bar we have γ∞ = 72mN/m.
The Tolman length has been estimated to be δT ≃
0.9Å from equilibrium measurements of the solvation
energy of spherical cavities.30 At the same conditions the
dynamic viscosity of SPC/E water has been found to
be η∞ = 6.42 · 10
−4Pa·s,31 ∼24% smaller than for real
water. In experiments in nanometer hydrophobic con-
finement and at interfaces however, the viscosity shows
deviations from the bulk value but remains comparable.26
We proceed by treating the viscosity η∞ as an adjustable
parameter together with its curvature correction coeffi-
cient δvis.
The MD simulations are carried out with the
DLPOLY2 package32 using an integration time step of
2fs. The simulation box is cubic and periodic in all three
dimensions with a length of L = (61.1 ± 0.2)Å in equi-
librium involving N = 6426 solvent molecules. Electro-
static interactions are calculated by the smooth-particle
mesh Ewald summation method. Lennard-Jones inter-
actions are cut-off and shifted at 9Å. Our investigated
systems are at first equilibrated in the NPT ensemble
with application of an external spherical potential of the
form βV (r) = [Å/(r −R′0)]
12 and all molecules removed
with r < R′0 since vapor can safely be neglected on these
scales. This stabilizes a well-defined spherical bubble of
radius R0 ≃ R
0 + 1Å. We define the cavity radius by
the radial location where the water density ρ(r) drops to
half of the bulk density ρ0/2. Thirty independent config-
urations in 20ps intervals are stored and serve as initial
configurations for the nonequilibrium runs. We employ a
Nosé-Hoover barostat and thermostat with a 0.2ps relax-
ation time to maintain the solvent at a pressure P and
a temperature T . Other choices of relaxation times in
the reasonable range between 0.1 and 0.5ps do not alter
our results. In the nonequilibrium simulations the con-
straining potential is switched off and the relaxation to
equilibrium is averaged over the thirty runs.
IV. RESULTS
system P/bar T/K cNaCl/M Q/e η∞/(10
−4Pa·s)
I 1 300 0 0 5.14
II 1 300 1.5 0 5.94
III 1 277 0 0 8.48
IV 2000 300 0 0 4.56
V 1000 300 0 0 4.72
VI 1 300 0 +2 5.14
TABLE I: Investigated system parameters: pressure P , tem-
perature T , and salt (NaCl) concentration c. In system VI
a fixed ion with charge Q = +2e is placed at the center of
the collapsing bubble. The viscosity η∞ is a fit-parameter in
systems I-V (see text).
We perform simulations of six different systems I-VI
whose features are summarized in Tab. I and differ in
thermodynamic parameters T and P (I, III, IV, and V)
but also inclusion of dispersed salt (II), and the influ-
ence of a charged particle in the bubble center (VI) are
considered. Note that the exact value of the crossover
length scale (however defined) can depend on the detailed
thermodynamic or solvent condition but remains close to
1nm.2
0 5 10 15 20 25 30
t=1ps
t=5ps
t=10ps
t=14ps
t=17ps
t=19ps
t=23ps
0 5 10 15 20
 τ(R)/
"10-90"- thickness
FIG. 1: Interface density profiles ρ(r)/ρ0 for system I are
plotted vs. the radial distance r from the bubble center for
different times t/ps=1,5,10,14,17,19,23. Symbols denote MD
simulation data and lines are fits using 2ρ(r)/ρ0 = erf{[r −
R(t)]/d}+1. The bubble radius R(t) is defined by the distance
at which the density is ρ0/2 (dotted line). The inset shows
the “10-90” thickness τ = 1.8124 d of the interface vs. R for
initial radii R0 = 19.83Å (pluses) and R0 = 25.6Å (crosses).
System I is at normal conditions (T=300K, P=1bar)
and consists of pure SPC/E water. Fig. 1 shows the
observed interface profiles in the nonequilibrium situa-
tion at different times t/ps=1, 5, 10, 14, 17, 19, and 23
starting from an initial radius R0 = 19.83Å. The liquid-
vapor interface stays relatively sharp in the process of
relaxation but broadens noticeably for smaller radii. At
t ≃ 23ps the system is completely relaxed to a homoge-
neous density distribution. The same time scale of bubble
collapse has been found in explicit water computer simu-
lations of dewetting in nanometer-sized paraffin plates,17
polymers,11 and atomistically resolved proteins.9,10
We find that the interface profiles can be fitted very
well with a functional form 2ρ(r)/ρ0 = erf{[r−R(t)]/d}+
1, where d is a measure of the interface thickness. The
interface fits are also shown in Fig. 1 together with the
MD data. The experimentally accessible “10-90” thick-
ness τ of an interface is the thickness over which the
density changes from 0.1ρ0 to 0.9ρ0 and is related to the
parameter d via τ = 1.8124 d. While experimental values
of τ for the planar water liquid-vapor interface vary be-
tween ∼ 4 and 8Å the measured values for SPC/E water
in the finite simulation systems are τ∞ =3 to 4Å.
29 We
find a strongly radius-dependent function τ(R) plotted
in the inset to Fig. 1 for initial radii R0 = 19.83Å and
R0 = 25.6Å. For R ≃ R0 the thickness increases during
the following 5ps from the equilibrium value τ ≃ 3Å to
FIG. 2: Time evolution of the cavity radius R(t) for parame-
ters as defined in systems I-VI. The solution of the modified
RP equation (4) (lines) is plotted vs. MD data (symbols).
The inset shows the solution of the modified RP equation in-
cluding inertia terms, cf. lhs of (1), (dashed lines) compared
to eq. (4) for system I with initial radii R0 = 19.83Å and
R0 = 10.0Å.
about τ ≃ 4.5− 5Å independent of R0. While the exact
equilibrium thickness at t = 0 depends on the particu-
lar choice of the confining potential V (r) (e.g., a softer
potential might lead to a broader initial interface) this
suggest that 4.5-5Å is the typical interface thickness
for a bubble of 1nm size. Regarding the slope of the
curve one might speculate that τ(R → ∞) saturates to
the thickness τ∞ of the measured planar interface for
R0 → ∞. For R . 10Å the thickness increases twofold
during the relaxation to equilibrium. This broadening
might be attributed to increased density fluctuations and
the structural change of interfacial water in the system
when crossing from large to small length scales which has
been shown to happen in equilibrium at ∼ 1nm.1,2
In Fig. 2 we plot the time evolution of the bubble ra-
dius R(t) for all investigated systems. Let us first focus
on the simulation data of system I (circles). As antic-
ipated the bubble radius decreases initially in a linear
fashion while for smaller radii (R(t) . 10Å) the velocity
steadily increases. From the best fit of eq. (4) we find
a viscosity η∞ = 5.14 · 10
−4Pa·s and its curvature cor-
rection coefficient δvis = 4.4Å. Although investigating a
confined system with large interfaces the viscosity value
differs only 20% from the SPC/E bulk value. Further-
more, from our macroscopic point of view the MD data
show that high curvature decreases the viscosity and the
latter has to be curvature-corrected with a (positive) co-
efficient larger than the Tolman length δT. If the surface
tension decreased in a stronger fashion with curvature
than viscosity the collapse velocity would drop in qual-
itative disagreement with the simulation. The overall
behavior of R(t) and the collapse velocity of about ∼
1Å/ps agrees very well with the recent MD data of Lugli
and Zerbetto, who simulated the collapse of a 1nm sized
bubble in SPC water.23
The inset to Fig. 2 shows the solution of eq. (4) in-
cluding inertial terms [left hand side of (1)] to check the
assumption of overdamped dynamics. While inertial ef-
fects are indeed small but not completely negligible for
an initial radius R0 = 19.83Å they basically vanish for
R0 = 10Å. Interestingly, the inertial effects are not visible
in the MD simulation data at all. We attribute this obser-
vation to the finite and periodic simulation box which is
known to suppress long-ranged inertial (hydrodynamic)
effects.33
In the following we assume δvis to be independent of
the other parameters and treat only η∞ as adjustable
variable. In system II we add 175 salt pairs of sodium
chloride (NaCl) into the aqueous solution resulting in a
concentration of c ≃1.5M. The ion-SPC/E interaction
parameters are those used by Bhatt et al.34 who mea-
sured a linear increase of surface tension with NaCl con-
centration in agreement with experimental data. While
this increment for c = 1.5M is about small 2-3%, the vis-
cosity has been measured experimentally to increase by
approximately 18% at 298.15K.35 Indeed by comparing
the simulation data to the theory we find a 16% larger
viscosity η∞ = 5.94 ·10
−4Pa·s. A slower collapse velocity
has been found also in the MD simulations of Lugli and
Zerbetto in concentrated LiCl and CsCL solutions when
compared to pure water.23
In system III we investigate the effect of lowering the
temperature by simulating at T = 277K.While only a 5%
increase of the water surface tension (SPC/E and real
water) is estimated from available data29 the viscosity
depends strongly on temperature: the relative increase
has been reported to be between 55 − 75% for SPC/E
water (85% for real water).36 Inspecting the MD data
and considering the surface tension increase we find in-
deed a large decrease in viscosity of 65% with a best-fit
η∞ = 8.48 · 10
−4Pa·s. Both systems, II and III, show
that solvent viscosity has a substantial influence on bub-
ble dynamics as quantitatively described by our simple
analytical approach. In systems IV and V we return to
T = 300K but increase the pressure P by a factor of
2000 and 1000, respectively. Best fits provide viscosities
which are around 10% smaller than at normal conditions
in agreement with the very weak pressure dependence of
the viscosity found in experiments37,38 at T=300K. The
major contribution to the faster dynamics comes explic-
itly from the pressure terms in eq. (4). Although mov-
ing away from liquid-vapor coexistence by increasing the
pressure up to 2000bar we assume (and verify hereby)
that the bubble interface tension can still be described
by γ∞.
In system VI we investigate the influence of a hy-
drophilic solute on the bubble interface motion in order
to make connection to cavitation close to molecular (pro-
tein) surfaces. As a simple measure we fix a divalent ion
at the center of the bubble so that we retain spherical
symmetry. The ion is modeled by a Lennard-Jones (LJ)
potential ULJ(r) = 4ǫ[(σ/r)
12 − (σ/r)6] with Q = +2e
point charges and uses the LJ parameters of the SPC/E
oxygen-oxygen interaction. As demonstrated recently
the LY equation can be modified to include dispersion
and electrostatic solute-solvent interactions explicitly,18
which extends (4) to
Ṙ = −
− ρ0ULJ(R) +
32πǫ0R4
The last term in (5) is the Born electrostatic energy den-
sity of a central charge Q in a spherical cavity with ra-
dius R with low dielectric vapor ǫv = 1 surrounded by
a high dielectric liquid (1/ǫl ≃ 0). The electric field
around the ionic charge and the dispersion attracts the
surrounding dipolar water what accelerates and eventu-
ally completely governs the bubble collapse below a ra-
dius R(t) . 13Å (t & 7ps) as also shown in Fig. 2.
The theoretical prediction (5) agrees very well without
any fitting using the viscosity from system I. We find
that the acceleration is mainly due to the electrostatic
attraction; the dispersion term plays just a minor role
while the excluded volume repulsion eventually deter-
mines the final (equilibrium) radius of the interface with
R(t = ∞) ≃ 2Å.
V. CONCLUSIONS
In conclusion, we have presented a simple analytical
and quantitative description of the interface motion of
a microscopic cavity by modifying the macroscopic RP
equation. Based on our MD data we find for the macro-
scopic description that analogous to the surface tension
the viscosity has to be corrected for curvature effects, a
prediction compelling to investigate further in detail and
probably related to the restructuring of interfacial water
for high curvatures (small R). The viscosity correction
accelerates collapse dynamics markedly below the equi-
librium crossover scale (∼1nm) in contrast to the pure
equilibrium picture where surface tension decreases what
slows down the collapse. Further, we find that the dy-
namics is curvature-driven due to the corrections to sur-
face tension and viscosity, not due surface tension as often
postulated.22 As a simple estimate, the interface velocity
is typically given by the ratio of surface tension and fluid
viscosity, v ≃ γ/η.
A comment has to be made regarding the recent work
of Lugli and Zerbetto on MD simulations of nanobubble
collapse in ionic solutions. While their MD data of the
collapse velocity for a 1nm bubble agree very well with
our results their interpretation in terms of the RP equa-
tion is different. They fit the ’violent regime’ solution
of the RP equation to the data [which is the solution of
only the inertial part, left hand side of (1)] and argue
that the violent regime still holds on the nm scale. As
demonstrated in this work, we arrive to a different con-
clusion: the collapse is friction dominated, the collapse
driving force is mainly capillary pressure, and we suggest
that the microscopic viscosity has to be curvature cor-
rected to explain the high curvature collapse behavior in
the MD simulations. The good agreement between our
modified RP equation and the MD data for different sol-
vent conditions, leading for instance to an altered solvent
surface tension or viscosity, support our view.
We finally note that extensions of the LY equation are
based on minimizing an appropriate free energy G(R) or
free energy functional18,24 so that we can write in a more
general form Ṙ ∼ [∂G(R)/∂R]/[η(R)R]. It is highly de-
sirable to generalize this simple dynamics further to ar-
bitrary geometries with which a wide field of potential
applications might open up, i.e. an efficient implicit mod-
eling of the water interface dynamics in the nonequilib-
rium process of hydrophobic nanoassembly, protein dock-
ing and folding, and nanofluidics.
Acknowledgements
J. D. thanks Lyderic Bocquet for pointing to the RP
equation, Bo Li (Applied Math, UCSD), Roland R.
Netz, Rudi Podgornik, and Dominik Horinek for stimu-
lating discussions, and the Deutsche Forschungsgemein-
schaft (DFG) for support within the Emmy-Noether-
Programme.
∗ e-mail address:jdzubiel@ph.tum.de
1 D. Chandler, Nature 437, 640 (2005).
2 S. Rajamani, T. M. Truskett, and S. Garde, PNAS 102,
9475 (2005).
3 A. Poynor, L. Hong, I. K. Robinson, S. Granick, Z. Zhang,
and P. A. Fenter, Phys. Rev. Lett. 97, 266101 (2006).
4 A. Wallquist and B. J. Berne, J. Phys. Chem. 99, 2893
(1995).
5 T. Koishi, Phys. Rev. Lett. 93, 185701 (2004).
6 N. Giovambattista, P. J. Rossky, and P. D. Debenedetti,
Phys. Rev. E 73, 041604 (2006).
7 O. Beckstein and M. S. P. Sansom, Proc. Nat. Acad. Sci.
100, 7063 (2003).
8 A. Anishkin and S. Sukharev, Biophys. J. 86, 2883 (2004).
9 R. Zhou, X. Huang, C. Margulis, and B. J. Berne, Science
305, 1605 (2004).
10 P. Liu, X. Huang, R. Zhou, and B. J. Berne, Nature 437,
159 (2005).
11 P. R. ten Wolde and D. Chandler, Proc. Natl. Acad. Sci.
99, 6539 (2002).
12 D. M. Huang and D. Chandler, PNAS 97, 8324 (2000).
13 A. Carambassis, L. C. Jonker, P. Attard, and M. W. Rut-
mailto:jdzubiel@ph.tum.de
land, Phys. Rev. Lett. 80, 5357 (1998).
14 R. Helmy, Y. Kazakevich, C. Ni, and A. Y. Fadeev, J. Am.
Chem. Soc. 127, 12446 (2005).
15 K. Jayaraman, K. Okamoto, S. J. Son, C. Luckett, A. H.
Gopalani, S. B. Lee, and D. S. English, J. Am. Chem. Soc.
127, 17385 (2005).
16 M. D. Collins, G. Hummer, M. L. Quillin, B. W. Matthews,
and S. M. Gruner, PNAS 102, 16668 (2005).
17 X. Huang, C. J. Margulis, and B. J. Berne, PNAS 100,
11953 (2003).
18 J. Dzubiella, J. M. J. Swanson, and J. A. McCammon,
Phys. Rev. Lett. 96, 087802 (2006).
19 J. Dzubiella, J. M. J. Swanson, and J. A. McCammon, J.
Chem. Phys. 124, 084905 (2006).
20 C. E. Brennen, Cavitation and bubble dynamics (Oxford
U. Press, New York, 1995).
21 M. S. Plesset and A. Prosperetti, Annu. Rev. Fluid. Mech.
25, 577 (1977).
22 H. Spohn, J. Stat. Phys. 71, 1081 (1993).
23 F. Lugli and F. Zerbetto, ChemPhysChem 8, 47 (2007).
24 L. Boruvka and A. W. Neumann, J. Chem. Phys. 66, 5464
(1977).
25 R. C. Tolman, J. Chem. Phys. 17, 333 (1949).
26 U. Raviv, S. Giasson, J. Frey, and J. Klein, J. Phys.: Con-
dens. Matt. 14, 9275 (2002).
27 I. W. Kuo and C. J. Mundy, Science 303, 658 (2004).
28 H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J.
Phys. Chem. 91, 6269 (1987).
29 J. Alejandre, D. J. Tildesley, and G. A. Chapela, J. Chem.
Phys. 102, 4574 (1995).
30 D. M. Huang and D. Chandler, J. Phys. Chem. B 106,
2047 (2002).
31 B. Hess, J. Chem. Phys. 116, 209 (2002).
32 W. Smith and T. R. Forester (1999), the DLPOLY 2 User
Manual.
33 B. Dünweg and K. Kremer, J. Chem. Phys. 99, 6983
(1993).
34 D. Bhatt, J. Newman, and C. J. Radke, J. Phys. Chem. B
108, 9077 (2004).
35 Z. Hai-Lang and H. Shi-Jun, J. Chem. Eng. Data 41, 516
(1996).
36 P. E. Smith and F. van Gunsteren, Chem. Phys. Lett. 215,
315 (1993).
37 K. E. Bett and J. B. Cappi, Nature 207, 620 (1965).
38 J. V. Sengers and J. T. R. Watson, J. Phys. Chem. Ref.
Data 15, 1291 (1986).

  An analytical description of the interface motion of a collapsing
nanometer-sized spherical cavity in water is presented by a modification of the
Rayleigh-Plesset equation in conjunction with explicit solvent molecular
dynamics simulations. Quantitative agreement is found between the two
approaches for the time-dependent cavity radius $R(t)$ at different solvent
conditions while in the continuum picture the solvent viscosity has to be
corrected for curvature effects. The typical magnitude of the interface or
collapse velocity is found to be given by the ratio of surface tension and
fluid viscosity, $v\simeq\gamma/\eta$, while the curvature correction
accelerates collapse dynamics on length scales below the equilibrium crossover
scales ($\sim$1nm). The study offers a starting point for an efficient implicit
modeling of water dynamics in aqueous nanoassembly and protein systems in
nonequilibrium.

Interface dynamics of microscopic cavities in water
Joachim Dzubiella1, ∗
Physics Department, Technical University Munich, 85748 Garching, Germany
(Dated: November 4, 2018)
An analytical description of the interface motion of a collapsing nanometer-sized spherical cavity
in water is presented by a modification of the Rayleigh-Plesset equation in conjunction with ex-
plicit solvent molecular dynamics simulations. Quantitative agreement is found between the two
approaches for the time-dependent cavity radius R(t) at different solvent conditions while in the con-
tinuum picture the solvent viscosity has to be corrected for curvature effects. The typical magnitude
of the interface or collapse velocity is found to be given by the ratio of surface tension and fluid vis-
cosity, v ≃ γ/η, while the curvature correction accelerates collapse dynamics on length scales below
the equilibrium crossover scales (∼1nm). The study offers a starting point for an efficient implicit
modeling of water dynamics in aqueous nanoassembly and protein systems in nonequilibrium.
I. INTRODUCTION
Hydrophobic hydration in equilibrium is a phe-
nomenon which exhibits qualitatively different behavior
at small and large length scales.1,2 While small solutes
(radii R .1nm) are accommodated by water with only
minor perturbations, larger solutes (R &1nm) induce ma-
jor rearrangements of water interfacial structure. As a
consequence the solvation free energy G(R) of small hy-
drophobic cavities scales with solute volume while for
larger cavities it grows with surface area (as a good
approximation near liquid-vapor coexistence) accompa-
nied by weak solvent dewetting at extended restrain-
ing hydrophobic surfaces.3 Furthermore, water, which
is close to the liquid-vapor transition at normal condi-
tions, can minimize interface area by locally evaporat-
ing and forming a ’nanobubble’ within hydrophobic con-
finement. Evidence of bubble formation in confined ge-
ometry has been given early by computer simulations
of smooth plate-like solutes,4 but more recently it has
been demonstrated in varying degrees in atomistically
resolved plate-like solutes,5,6 hydrophobic tubes and ion
channels,7,8 and in the collapse of proteins,9,10 suggesting
that it plays a key role in the stabilization and folding
dynamics of certain classes of biomolecules.11,12 Experi-
mental evidence of nanobubbles in strong confinement (in
contrast to bubbles at a single planar surface3) has been
given for instance in studies of water between hydropho-
bic surfaces,13 in zeolites and silica nanotubes,14,15 and
on a subnanometer scale in nonpolar protein cavities.16
The dewetting induced change in solvation energy is
typically estimated using simple macroscopic arguments
as known from capillarity theory, e.g. by describing in-
terfaces with Laplace-Young (LY) type of equations.14,17
Recently an extension of the LY equation has become
available which extrapolates to microscopic scales by in-
cluding a curvature correction to the interface tension
and considering atomistic dispersion and electrostatic
potentials of the solvated solute explicitly.18 Although
those macroscopic considerations (e.g,. the concept of
surface tension) are supposed to break down on atom-
istic scales they show surprisingly good results for the
solvation energy of microscopic solutes, e.g. alkanes and
noble gases, and quantitatively account for dewetting
effects in nanometer-sized hydrophobic confinement.19
While we conclude that the equilibrium location of the
solute-solvent interface seems to be well described by
those techniques, nothing is known about the interface
dynamics of evolution and relaxation. In this study we
address two fundamental questions: First, what are the
equations which govern the interface motion on atomistic
(∼1nm) scales? Secondly, does the dynamics exhibit any
signatures of the length scale crossover found in equilib-
rium?
On macroscopic scales the collapse dynamics of a (va-
por or gas) bubble is related to the well-known phe-
nomenon of sonoluminescence.20 The governing equa-
tions can be derived from Navier-Stokes and capillarity
theory and are expressed by the Rayleigh-Plesset (RP)
equation.21 We will show that the RP equation simpli-
fies in the limit of microscopic cavities and can be ex-
tended to give a quantitative description of cavity inter-
face dynamics on nanometer length scales. We find a
qualitatively different dynamics than the typical “mean-
curvature flow” description of moving interfaces,22 in par-
ticular a typical magnitude of interface or collapse veloc-
ity given by the ratio of surface tension and fluid viscosity,
v ≃ γ/η. Our study is restricted to the generic case of
the collapse of a spherical cavity and is complemented
by explicit solvent molecular dynamics (MD) computer
simulations. We note here that recently, Lugli and Zer-
betto studied nanobubble collapse in ionic solutions by
MD simulations on similar length scales.23 While their
MD data compares favorably with our results their in-
terpretation and conclusions in terms of the RP equa-
tion are different. We will resume this discussion in the
conclusion section.
In this study we show that a simple analytical approach
quantitatively describes microscopic cavity collapse for
a variety of different solvent situations while the sim-
ulations suggest that the solvent viscosity needs to be
corrected for curvature effects. Our study might offer a
simple starting point for an efficient implicit modeling
of water dynamics in aqueous nanoassembly and protein
http://arxiv.org/abs/0704.0239v1
systems in nonequilibrium.
II. THEORY
The Rayleigh-Plesset equation for the time evolution
of a macroscopic vapor bubble with radius R(t) can be
written as21
= ∆P + 4η
, (1)
where ρm is the solvent mass density, ∆P = P − Pv the
difference in liquid and vapor pressures, η the dynamic
viscosity, and γ the liquid-vapor interface tension. While
for macroscopic bubble radii the inertial terms (left hand
side) control the dynamics, for decreasing radii the fric-
tional and pressure terms (right hand side) grow in rel-
ative magnitude and eventually dominate, so that com-
pletely overdamped dynamics can be assumed on atom-
istic scales:
Ṙ ≃ −
. (2)
A rough estimate for the threshold radius Rt below which
friction dominates is given when the Reynolds number
R = vRρm/η becomes unity and viscous and inertial
forces are balanced. With a typical initial interface ve-
locity of the order of v ∼ γ/η [from R̈(0) = 0 in eq. (1)]
we obtain
Rt = η
2/(ρmγ) (3)
which is ≃ 10nm for water at normal conditions. Note
that this threshold value can deviate considerably for a
fluid different than water and that the viscosity typically
has a strong temperature (T ) dependence which implies
that Rt can change significantly with T .
In equilibrium (Ṙ = 0) the remaining expression in
eq. (2) is the (spherical) LY equation ∆P + 2γ/R = 0.
Thus eq. (2) describes a linear relationship between cap-
illary pressure and interface velocity where R/(4η) plays
the role of an interface mobility (inverse friction).22 In-
terestingly, the mobility is linear in bubble radius which
leads to a constant velocity driven by surface tension in-
dependent of radius (assuming P ≃ 0); this has to be con-
trasted to the typically used capillary dynamics which is
proportional to the local mean curvature ∝ 1/R.22
Generalizations of the LY equation to small scales are
available by adding a Gaussian curvature term (∼ 1/R2)
as shown by Boruvka and Neumann24; that has been
demonstrated to be equivalent to a first order curvature
correction in surface tension, i.e. γ(R) = γ∞(1−δT/R),
where δT is the Tolman length
25 and γ∞ the liquid-vapor
surface tension for a planar interface (R = ∞). The
Tolman length has a magnitude which is usually of the
order of the size of a solvent molecule. Furthermore, it
has been observed experimentally that the viscosity of
strongly confined water can depend on the particular na-
ture of the surface/interface.26 We conclude that in gen-
eral one has to anticipate that - analogous to the surface
tension - the effective interface viscosity obeys a curva-
ture correction in the limit of small cavities due to water
restructuring in the first solvent layers at the hydropho-
bic interface. In the following we make the simple first
order assumption that the correction enters eq. (2) also
linear in curvature (∼ 1/R) yielding
Ṙ = −
∆PR+∆Pδvis + 2γ∞ +
, (4)
where the constant δvis is the coefficient for the first order
curvature correction in viscosity and η∞ the macroscopic
bulk viscosity. Additionally, we define δ = δvis − δT and
second order terms in curvature are neglected. We note
that the choice of the 1/R-scaling of the viscosity curva-
ture correction has no direct physical justification and is
arbitrary. We think however, that a curvature correction
based on an expansion in orders of mean curvature is the
simplest and most natural way for such a choice.
In water at normal conditions the pressure terms in
(4) are negligible so that for large radii (R ≫ δ) the in-
terface velocity is constant and R(t) = R0 − γ∞/(2η∞)t.
This leads to a collapse velocity of about v ≃0.4Å/ps
(40m/s) which is 6% of the thermal velocity of water
vth =
3kBT/m showing that dissipative heating of the
system is relatively weak on these scales. A rough esti-
mate for the dissipation rate can be made by the released
interfacial energy dG(R, t)/dt ≃ d(4πR(t)2γ∞)/dt =
−4πγ2
R(t)/η∞ yielding for instance dG(R, t = 0)/dt ≃
−35kBT/ps for a bubble with R0 =2nm. At small
radii (R ≃ δ) the solution of (4) goes as R(t) ∼
const− (δγ∞/η∞)t decreasing or increasing the ve-
locity depending on the sign of δ = δvis − δT, i.e. the
acceleration depends on the particular sign and mag-
nitude of the curvature corrections to surface tension
and viscosity. For large pressures and radii the first
term dominates which gives rise to an exponential de-
cay R(t) ∼ exp[−∆P/(4η∞)t]. While extending to small
scales we have assumed that the time scale of internal in-
terface dynamics, i.e. hydrogen bond rearrangements,27
is much faster than the one of bubble collapse.
III. MD SIMULATION
In order to quantify our analytical predictions we
complement the theory by MD simulations using ex-
plicit SPC/E water.28 The liquid-vapor surface tension
of SPC/E water has been measured and agrees with the
experimental value for a wide range of temperatures.29
For T = 300K and P = 1bar we have γ∞ = 72mN/m.
The Tolman length has been estimated to be δT ≃
0.9Å from equilibrium measurements of the solvation
energy of spherical cavities.30 At the same conditions the
dynamic viscosity of SPC/E water has been found to
be η∞ = 6.42 · 10
−4Pa·s,31 ∼24% smaller than for real
water. In experiments in nanometer hydrophobic con-
finement and at interfaces however, the viscosity shows
deviations from the bulk value but remains comparable.26
We proceed by treating the viscosity η∞ as an adjustable
parameter together with its curvature correction coeffi-
cient δvis.
The MD simulations are carried out with the
DLPOLY2 package32 using an integration time step of
2fs. The simulation box is cubic and periodic in all three
dimensions with a length of L = (61.1 ± 0.2)Å in equi-
librium involving N = 6426 solvent molecules. Electro-
static interactions are calculated by the smooth-particle
mesh Ewald summation method. Lennard-Jones inter-
actions are cut-off and shifted at 9Å. Our investigated
systems are at first equilibrated in the NPT ensemble
with application of an external spherical potential of the
form βV (r) = [Å/(r −R′0)]
12 and all molecules removed
with r < R′0 since vapor can safely be neglected on these
scales. This stabilizes a well-defined spherical bubble of
radius R0 ≃ R
0 + 1Å. We define the cavity radius by
the radial location where the water density ρ(r) drops to
half of the bulk density ρ0/2. Thirty independent config-
urations in 20ps intervals are stored and serve as initial
configurations for the nonequilibrium runs. We employ a
Nosé-Hoover barostat and thermostat with a 0.2ps relax-
ation time to maintain the solvent at a pressure P and
a temperature T . Other choices of relaxation times in
the reasonable range between 0.1 and 0.5ps do not alter
our results. In the nonequilibrium simulations the con-
straining potential is switched off and the relaxation to
equilibrium is averaged over the thirty runs.
IV. RESULTS
system P/bar T/K cNaCl/M Q/e η∞/(10
−4Pa·s)
I 1 300 0 0 5.14
II 1 300 1.5 0 5.94
III 1 277 0 0 8.48
IV 2000 300 0 0 4.56
V 1000 300 0 0 4.72
VI 1 300 0 +2 5.14
TABLE I: Investigated system parameters: pressure P , tem-
perature T , and salt (NaCl) concentration c. In system VI
a fixed ion with charge Q = +2e is placed at the center of
the collapsing bubble. The viscosity η∞ is a fit-parameter in
systems I-V (see text).
We perform simulations of six different systems I-VI
whose features are summarized in Tab. I and differ in
thermodynamic parameters T and P (I, III, IV, and V)
but also inclusion of dispersed salt (II), and the influ-
ence of a charged particle in the bubble center (VI) are
considered. Note that the exact value of the crossover
length scale (however defined) can depend on the detailed
thermodynamic or solvent condition but remains close to
1nm.2
0 5 10 15 20 25 30
t=1ps
t=5ps
t=10ps
t=14ps
t=17ps
t=19ps
t=23ps
0 5 10 15 20
 τ(R)/
"10-90"- thickness
FIG. 1: Interface density profiles ρ(r)/ρ0 for system I are
plotted vs. the radial distance r from the bubble center for
different times t/ps=1,5,10,14,17,19,23. Symbols denote MD
simulation data and lines are fits using 2ρ(r)/ρ0 = erf{[r −
R(t)]/d}+1. The bubble radius R(t) is defined by the distance
at which the density is ρ0/2 (dotted line). The inset shows
the “10-90” thickness τ = 1.8124 d of the interface vs. R for
initial radii R0 = 19.83Å (pluses) and R0 = 25.6Å (crosses).
System I is at normal conditions (T=300K, P=1bar)
and consists of pure SPC/E water. Fig. 1 shows the
observed interface profiles in the nonequilibrium situa-
tion at different times t/ps=1, 5, 10, 14, 17, 19, and 23
starting from an initial radius R0 = 19.83Å. The liquid-
vapor interface stays relatively sharp in the process of
relaxation but broadens noticeably for smaller radii. At
t ≃ 23ps the system is completely relaxed to a homoge-
neous density distribution. The same time scale of bubble
collapse has been found in explicit water computer simu-
lations of dewetting in nanometer-sized paraffin plates,17
polymers,11 and atomistically resolved proteins.9,10
We find that the interface profiles can be fitted very
well with a functional form 2ρ(r)/ρ0 = erf{[r−R(t)]/d}+
1, where d is a measure of the interface thickness. The
interface fits are also shown in Fig. 1 together with the
MD data. The experimentally accessible “10-90” thick-
ness τ of an interface is the thickness over which the
density changes from 0.1ρ0 to 0.9ρ0 and is related to the
parameter d via τ = 1.8124 d. While experimental values
of τ for the planar water liquid-vapor interface vary be-
tween ∼ 4 and 8Å the measured values for SPC/E water
in the finite simulation systems are τ∞ =3 to 4Å.
29 We
find a strongly radius-dependent function τ(R) plotted
in the inset to Fig. 1 for initial radii R0 = 19.83Å and
R0 = 25.6Å. For R ≃ R0 the thickness increases during
the following 5ps from the equilibrium value τ ≃ 3Å to
FIG. 2: Time evolution of the cavity radius R(t) for parame-
ters as defined in systems I-VI. The solution of the modified
RP equation (4) (lines) is plotted vs. MD data (symbols).
The inset shows the solution of the modified RP equation in-
cluding inertia terms, cf. lhs of (1), (dashed lines) compared
to eq. (4) for system I with initial radii R0 = 19.83Å and
R0 = 10.0Å.
about τ ≃ 4.5− 5Å independent of R0. While the exact
equilibrium thickness at t = 0 depends on the particu-
lar choice of the confining potential V (r) (e.g., a softer
potential might lead to a broader initial interface) this
suggest that 4.5-5Å is the typical interface thickness
for a bubble of 1nm size. Regarding the slope of the
curve one might speculate that τ(R → ∞) saturates to
the thickness τ∞ of the measured planar interface for
R0 → ∞. For R . 10Å the thickness increases twofold
during the relaxation to equilibrium. This broadening
might be attributed to increased density fluctuations and
the structural change of interfacial water in the system
when crossing from large to small length scales which has
been shown to happen in equilibrium at ∼ 1nm.1,2
In Fig. 2 we plot the time evolution of the bubble ra-
dius R(t) for all investigated systems. Let us first focus
on the simulation data of system I (circles). As antic-
ipated the bubble radius decreases initially in a linear
fashion while for smaller radii (R(t) . 10Å) the velocity
steadily increases. From the best fit of eq. (4) we find
a viscosity η∞ = 5.14 · 10
−4Pa·s and its curvature cor-
rection coefficient δvis = 4.4Å. Although investigating a
confined system with large interfaces the viscosity value
differs only 20% from the SPC/E bulk value. Further-
more, from our macroscopic point of view the MD data
show that high curvature decreases the viscosity and the
latter has to be curvature-corrected with a (positive) co-
efficient larger than the Tolman length δT. If the surface
tension decreased in a stronger fashion with curvature
than viscosity the collapse velocity would drop in qual-
itative disagreement with the simulation. The overall
behavior of R(t) and the collapse velocity of about ∼
1Å/ps agrees very well with the recent MD data of Lugli
and Zerbetto, who simulated the collapse of a 1nm sized
bubble in SPC water.23
The inset to Fig. 2 shows the solution of eq. (4) in-
cluding inertial terms [left hand side of (1)] to check the
assumption of overdamped dynamics. While inertial ef-
fects are indeed small but not completely negligible for
an initial radius R0 = 19.83Å they basically vanish for
R0 = 10Å. Interestingly, the inertial effects are not visible
in the MD simulation data at all. We attribute this obser-
vation to the finite and periodic simulation box which is
known to suppress long-ranged inertial (hydrodynamic)
effects.33
In the following we assume δvis to be independent of
the other parameters and treat only η∞ as adjustable
variable. In system II we add 175 salt pairs of sodium
chloride (NaCl) into the aqueous solution resulting in a
concentration of c ≃1.5M. The ion-SPC/E interaction
parameters are those used by Bhatt et al.34 who mea-
sured a linear increase of surface tension with NaCl con-
centration in agreement with experimental data. While
this increment for c = 1.5M is about small 2-3%, the vis-
cosity has been measured experimentally to increase by
approximately 18% at 298.15K.35 Indeed by comparing
the simulation data to the theory we find a 16% larger
viscosity η∞ = 5.94 ·10
−4Pa·s. A slower collapse velocity
has been found also in the MD simulations of Lugli and
Zerbetto in concentrated LiCl and CsCL solutions when
compared to pure water.23
In system III we investigate the effect of lowering the
temperature by simulating at T = 277K.While only a 5%
increase of the water surface tension (SPC/E and real
water) is estimated from available data29 the viscosity
depends strongly on temperature: the relative increase
has been reported to be between 55 − 75% for SPC/E
water (85% for real water).36 Inspecting the MD data
and considering the surface tension increase we find in-
deed a large decrease in viscosity of 65% with a best-fit
η∞ = 8.48 · 10
−4Pa·s. Both systems, II and III, show
that solvent viscosity has a substantial influence on bub-
ble dynamics as quantitatively described by our simple
analytical approach. In systems IV and V we return to
T = 300K but increase the pressure P by a factor of
2000 and 1000, respectively. Best fits provide viscosities
which are around 10% smaller than at normal conditions
in agreement with the very weak pressure dependence of
the viscosity found in experiments37,38 at T=300K. The
major contribution to the faster dynamics comes explic-
itly from the pressure terms in eq. (4). Although mov-
ing away from liquid-vapor coexistence by increasing the
pressure up to 2000bar we assume (and verify hereby)
that the bubble interface tension can still be described
by γ∞.
In system VI we investigate the influence of a hy-
drophilic solute on the bubble interface motion in order
to make connection to cavitation close to molecular (pro-
tein) surfaces. As a simple measure we fix a divalent ion
at the center of the bubble so that we retain spherical
symmetry. The ion is modeled by a Lennard-Jones (LJ)
potential ULJ(r) = 4ǫ[(σ/r)
12 − (σ/r)6] with Q = +2e
point charges and uses the LJ parameters of the SPC/E
oxygen-oxygen interaction. As demonstrated recently
the LY equation can be modified to include dispersion
and electrostatic solute-solvent interactions explicitly,18
which extends (4) to
Ṙ = −
− ρ0ULJ(R) +
32πǫ0R4
The last term in (5) is the Born electrostatic energy den-
sity of a central charge Q in a spherical cavity with ra-
dius R with low dielectric vapor ǫv = 1 surrounded by
a high dielectric liquid (1/ǫl ≃ 0). The electric field
around the ionic charge and the dispersion attracts the
surrounding dipolar water what accelerates and eventu-
ally completely governs the bubble collapse below a ra-
dius R(t) . 13Å (t & 7ps) as also shown in Fig. 2.
The theoretical prediction (5) agrees very well without
any fitting using the viscosity from system I. We find
that the acceleration is mainly due to the electrostatic
attraction; the dispersion term plays just a minor role
while the excluded volume repulsion eventually deter-
mines the final (equilibrium) radius of the interface with
R(t = ∞) ≃ 2Å.
V. CONCLUSIONS
In conclusion, we have presented a simple analytical
and quantitative description of the interface motion of
a microscopic cavity by modifying the macroscopic RP
equation. Based on our MD data we find for the macro-
scopic description that analogous to the surface tension
the viscosity has to be corrected for curvature effects, a
prediction compelling to investigate further in detail and
probably related to the restructuring of interfacial water
for high curvatures (small R). The viscosity correction
accelerates collapse dynamics markedly below the equi-
librium crossover scale (∼1nm) in contrast to the pure
equilibrium picture where surface tension decreases what
slows down the collapse. Further, we find that the dy-
namics is curvature-driven due to the corrections to sur-
face tension and viscosity, not due surface tension as often
postulated.22 As a simple estimate, the interface velocity
is typically given by the ratio of surface tension and fluid
viscosity, v ≃ γ/η.
A comment has to be made regarding the recent work
of Lugli and Zerbetto on MD simulations of nanobubble
collapse in ionic solutions. While their MD data of the
collapse velocity for a 1nm bubble agree very well with
our results their interpretation in terms of the RP equa-
tion is different. They fit the ’violent regime’ solution
of the RP equation to the data [which is the solution of
only the inertial part, left hand side of (1)] and argue
that the violent regime still holds on the nm scale. As
demonstrated in this work, we arrive to a different con-
clusion: the collapse is friction dominated, the collapse
driving force is mainly capillary pressure, and we suggest
that the microscopic viscosity has to be curvature cor-
rected to explain the high curvature collapse behavior in
the MD simulations. The good agreement between our
modified RP equation and the MD data for different sol-
vent conditions, leading for instance to an altered solvent
surface tension or viscosity, support our view.
We finally note that extensions of the LY equation are
based on minimizing an appropriate free energy G(R) or
free energy functional18,24 so that we can write in a more
general form Ṙ ∼ [∂G(R)/∂R]/[η(R)R]. It is highly de-
sirable to generalize this simple dynamics further to ar-
bitrary geometries with which a wide field of potential
applications might open up, i.e. an efficient implicit mod-
eling of the water interface dynamics in the nonequilib-
rium process of hydrophobic nanoassembly, protein dock-
ing and folding, and nanofluidics.
Acknowledgements
J. D. thanks Lyderic Bocquet for pointing to the RP
equation, Bo Li (Applied Math, UCSD), Roland R.
Netz, Rudi Podgornik, and Dominik Horinek for stimu-
lating discussions, and the Deutsche Forschungsgemein-
schaft (DFG) for support within the Emmy-Noether-
Programme.
∗ e-mail address:jdzubiel@ph.tum.de
1 D. Chandler, Nature 437, 640 (2005).
2 S. Rajamani, T. M. Truskett, and S. Garde, PNAS 102,
9475 (2005).
3 A. Poynor, L. Hong, I. K. Robinson, S. Granick, Z. Zhang,
and P. A. Fenter, Phys. Rev. Lett. 97, 266101 (2006).
4 A. Wallquist and B. J. Berne, J. Phys. Chem. 99, 2893
(1995).
5 T. Koishi, Phys. Rev. Lett. 93, 185701 (2004).
6 N. Giovambattista, P. J. Rossky, and P. D. Debenedetti,
Phys. Rev. E 73, 041604 (2006).
7 O. Beckstein and M. S. P. Sansom, Proc. Nat. Acad. Sci.
100, 7063 (2003).
8 A. Anishkin and S. Sukharev, Biophys. J. 86, 2883 (2004).
9 R. Zhou, X. Huang, C. Margulis, and B. J. Berne, Science
305, 1605 (2004).
10 P. Liu, X. Huang, R. Zhou, and B. J. Berne, Nature 437,
159 (2005).
11 P. R. ten Wolde and D. Chandler, Proc. Natl. Acad. Sci.
99, 6539 (2002).
12 D. M. Huang and D. Chandler, PNAS 97, 8324 (2000).
13 A. Carambassis, L. C. Jonker, P. Attard, and M. W. Rut-
mailto:jdzubiel@ph.tum.de
land, Phys. Rev. Lett. 80, 5357 (1998).
14 R. Helmy, Y. Kazakevich, C. Ni, and A. Y. Fadeev, J. Am.
Chem. Soc. 127, 12446 (2005).
15 K. Jayaraman, K. Okamoto, S. J. Son, C. Luckett, A. H.
Gopalani, S. B. Lee, and D. S. English, J. Am. Chem. Soc.
127, 17385 (2005).
16 M. D. Collins, G. Hummer, M. L. Quillin, B. W. Matthews,
and S. M. Gruner, PNAS 102, 16668 (2005).
17 X. Huang, C. J. Margulis, and B. J. Berne, PNAS 100,
11953 (2003).
18 J. Dzubiella, J. M. J. Swanson, and J. A. McCammon,
Phys. Rev. Lett. 96, 087802 (2006).
19 J. Dzubiella, J. M. J. Swanson, and J. A. McCammon, J.
Chem. Phys. 124, 084905 (2006).
20 C. E. Brennen, Cavitation and bubble dynamics (Oxford
U. Press, New York, 1995).
21 M. S. Plesset and A. Prosperetti, Annu. Rev. Fluid. Mech.
25, 577 (1977).
22 H. Spohn, J. Stat. Phys. 71, 1081 (1993).
23 F. Lugli and F. Zerbetto, ChemPhysChem 8, 47 (2007).
24 L. Boruvka and A. W. Neumann, J. Chem. Phys. 66, 5464
(1977).
25 R. C. Tolman, J. Chem. Phys. 17, 333 (1949).
26 U. Raviv, S. Giasson, J. Frey, and J. Klein, J. Phys.: Con-
dens. Matt. 14, 9275 (2002).
27 I. W. Kuo and C. J. Mundy, Science 303, 658 (2004).
28 H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J.
Phys. Chem. 91, 6269 (1987).
29 J. Alejandre, D. J. Tildesley, and G. A. Chapela, J. Chem.
Phys. 102, 4574 (1995).
30 D. M. Huang and D. Chandler, J. Phys. Chem. B 106,
2047 (2002).
31 B. Hess, J. Chem. Phys. 116, 209 (2002).
32 W. Smith and T. R. Forester (1999), the DLPOLY 2 User
Manual.
33 B. Dünweg and K. Kremer, J. Chem. Phys. 99, 6983
(1993).
34 D. Bhatt, J. Newman, and C. J. Radke, J. Phys. Chem. B
108, 9077 (2004).
35 Z. Hai-Lang and H. Shi-Jun, J. Chem. Eng. Data 41, 516
(1996).
36 P. E. Smith and F. van Gunsteren, Chem. Phys. Lett. 215,
315 (1993).
37 K. E. Bett and J. B. Cappi, Nature 207, 620 (1965).
38 J. V. Sengers and J. T. R. Watson, J. Phys. Chem. Ref.
Data 15, 1291 (1986).
