# Molecular mechanism for cavitation in water under tension

Georg Menzl, Miguel A. Gonzalez, Philipp Geiger, Frederic Caupin, Jose L. F. Abascal, Chantal Valeriani, Christoph Dellago

► To cite this version:

Georg Menzl, Miguel A. Gonzalez, Philipp Geiger, Frederic Caupin, Jose L. F. Abascal, et al. Molec- ular mechanism for cavitation in water under tension. Proceedings of the National Academy of Sciences of the United States of America, 2016, 113, pp.13582-13587. 10.1073/pnas.1608421113 . hal-02290213

HAL Id: hal-02290213

https://univ-lyon1.hal.science/hal-02290213v1

Submitted on 31 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Molecular mechanism for cavitation in water under tension

Georg Menzl<sup>a,b</sup>, Miguel A. Gonzalez<sup>c</sup>, Philipp Geiger<sup>a,b</sup>, Frédéric Caupin<sup>d</sup>, José L. F. Abascal<sup>e</sup>, Chantal Valeriani<sup>e,f</sup>, and Christoph Dellago<sup>a,b,1</sup>

Edited by Daan Frenkel, University of Cambridge, Cambridge, United Kingdom, and approved September 23, 2016 (received for review May 25, 2016)

Despite its relevance in biology and engineering, the molecular mechanism driving cavitation in water remains unknown. Using computer simulations, we investigate the structure and dynamics of vapor bubbles emerging from metastable water at negative pressures. We find that in the early stages of cavitation, bubbles are irregularly shaped and become more spherical as they grow. Nevertheless, the free energy of bubble formation can be perfectly reproduced in the framework of classical nucleation theory (CNT) if the curvature dependence of the surface tension is taken into account. Comparison of the observed bubble dynamics to the predictions of the macroscopic Rayleigh-Plesset (RP) equation, augmented with thermal fluctuations, demonstrates that the growth of nanoscale bubbles is governed by viscous forces. Combining the dynamical prefactor determined from the RP equation with CNT based on the Kramers formalism yields an analytical expression for the cavitation rate that reproduces the simulation results very well over a wide range of pressures. Furthermore, our theoretical predictions are in excellent agreement with cavitation rates obtained from inclusion experiments. This suggests that homogeneous nucleation is observed in inclusions, whereas only heterogeneous nucleation on impurities or defects occurs in other experiments.

cavitation | water | negative pressure | bubble nucleation | liquid-vapor transition

Due to its pronounced cohesion, water remains stable under tension for long times. Experimentally, strongly negative pressures exceeding –120 MPa (1–6) can be sustained before the system decays into the vapor phase via cavitation, i.e., bubble nucleation. Recently, cavitation in water under tension has drawn research interest due to its importance in biological processes, like water transport in natural (7–10) and synthetic (11, 12) trees, spore propagation of ferns (13), and poration of cell membranes (14, 15). Furthermore, cavitation in water appears to be the driving force behind the sonocrystallization of ice (16, 17), and preventing its occurrence remains a challenge in turbine and propeller design (18). Studying the onset of cavitation has also proven to be a valuable tool to locate the line of density maxima in metastable water (4), which contributes to the ongoing effort of explaining the origin of water's anomalies (6, 19). Interest in the topic is magnified by the startling discrepancy arising when cavitation in water is investigated using different experimental methods. Although agreement between different methods is excellent in the high-temperature regime, where the liquid is unable to sustain large tension, a significantly higher degree of metastability is reached when studying cavitation in inclusions along an isochoric path (1–5) compared with other techniques (20, 21) at low temperatures (22).

Due to the short time scale on which the transition takes place and the small volume of the critical bubble at experimentally feasible conditions, direct observation of cavitation at the microscopic level remains elusive. However, cavitation rates are directly accessible in experiment and some microscopic insight into the cavitation transition can be obtained from these data by

means of the nucleation theorem (23), which relates the variation in the height of the free energy barrier separating the metastable liquid from the vapor phase upon change of external parameters to properties of the critical bubble (4, 21). The microscopic information that can be inferred is limited, and because not all quantities entering the nucleation theorem are known, ad hoc assumptions have to be introduced. For state points where cavitation is a rare event, classical nucleation theory (CNT) can be invoked to provide a qualitative understanding of the transition (24). However, although CNT provides a physically meaningful and appealingly simple picture of nucleation processes, the estimates for the nucleation rates obtained from CNT are known to differ substantially (up to many orders of magnitude) from those measured in experiments (22, 25, 26).

Computer simulations are a natural choice to investigate cavitation in water with molecular resolution on the time scales governing the emergence of microscopic bubbles in the liquid. Although cavitation in simple liquids has been studied extensively using computer simulations (27–33), simulation studies of cavitation in water were focused on methodological aspects (34–36) or performed at state points in vicinity of the vapor–liquid spinodal (37, 38). In this work, we apply a combination of several complementary computer simulation methods to identify the molecular mechanism of cavitation. A statistical committor analysis carried out on reactive trajectories reveals that the volume of the largest bubble in the system constitutes a good reaction coordinate

## Significance

Cavitation, the formation of vapor-filled bubbles in a liquid at low pressures, is a powerful phenomenon with important consequences in nature and technology. For instance, cavitation bubbles may interrupt water flow in plants under dry conditions or severely damage the metal surfaces of machines such as pumps and propellers. Using molecular simulations, we have studied cavitation in water at strongly negative pressures and have revealed its molecular mechanism. We find that bubble growth is governed by the viscosity of the liquid. Although small bubbles are shaped irregularly, classical nucleation theory accurately describes the free energy barrier that impedes rapid bubble formation. Our simulations indicate that water can withstand negative pressures exceeding –120 MPa in agreement with recent experiments.

Author contributions: G.M., F.C., C.V., and C.D. designed research; G.M., P.G., and C.D. performed research; G.M., M.A.G., F.C., J.L.F.A., C.V., and C.D. analyzed data; and G.M., M.A.G., P.G., F.C., J.L.F.A., C.V., and C.D. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

See Commentary on page 13549.

<footnote><sup>a</sup>Faculty of Physics, University of Vienna, 1090 Vienna, Austria; <sup>b</sup>Center for Computational Materials Science, University of Vienna, 1090 Vienna, Austria;</footnote>

<footnote><sup>c</sup>Department of Chemistry, Imperial College London, London SW7 2AZ, United Kingdom; <sup>d</sup>Institut Lumière Matière, UMR5306 Université Claude Bernard Lyon 1, CNRS, Université de Lyon, Institut Universitaire de France, 69622 Villeurbanne, France; <sup>e</sup>Departamento de Química Física, Facultad de Ciencias Químicas, Universidad Complutense de Madrid, 28040 Madrid, Spain; and <sup>f</sup>Departamento de Física Aplicada I, Facultad de Ciencias Física, Universidad Complutense de Madrid, 28040 Madrid, Spain</footnote>

<footnote n="1">¹To whom correspondence should be addressed. Email: christoph.dellago@univie.ac.at.</footnote>

<footnote>This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1608421113/-/DCSupplemental.</footnote>
for bubble nucleation. We compute the dynamics of nanoscale bubbles along this reaction coordinate and demonstrate that the pressure dependence of the bubble diffusivity can be reproduced by Rayleigh-Plesset (RP) theory generalized to include thermal fluctuations, thereby elucidating the crucial influence of viscous damping on bubble growth. Based on the Kramers formalism and the RP equation we obtain an analytical expression for the nucleation rate that yields excellent agreement with numerical results obtained for a wide range of pressures with a method akin to the Bennett-Chandler approach for the computation of reaction rate constants. The obtained rates are validated for selected points by comparison with estimates from transition interface sampling and support estimates obtained from inclusion experiments. To augment the microscopic picture of cavitation, we characterize the morphology of bubbles in water under tension and analyze the bubble surface in terms of its hydrogen bonding structure.

## Classical Nucleation Theory

Our investigations are guided by CNT, which posits that the decay of the metastable liquid under tension proceeds via the formation of a small vapor bubble, whose growth is initially opposed by a free energy barrier. According to Kramers theory (39, 40), the escape rate k from a well over a high barrier for a system moving diffusively in a potential U(q) along a coordinate q is given by $k = [(\int_0^\infty \exp[-\beta U(q)] dq)(\int_0^\infty \exp[\beta U(q)]/D(q) dq)]^{-1}$. Here the symbols ∪ and ∩ indicate that the integration is carried out over the well and the barrier, respectively, and D(q) is the diffusion coefficient. To describe bubble nucleation, we use the volume v of the largest bubble in the system as the order parameter [committor calculations (41) indicate that v is indeed a good reaction coordinate; Fig. S1], and we replace the potential energy by the potential of mean force $-k_B T \ln|v_0 P(v)|$, where k_B is the Boltzmann constant, P(v) is the probability density that the largest bubble is of size v, and v_0 is an arbitrary constant volume. Assuming that the diffusion coefficient does not change appreciably on the top of the barrier and approximating the barrier to second order, one obtains the nucleation rate (number of nucleation events per unit time and unit volume)

$$J = \frac{\omega D(v^*) P(v^*)}{\sqrt{2\pi k_B T} V}, \quad [1]$$

where v^* is the critical bubble volume, V is the total volume of the system, and ω is related to the barrier curvature κ by ω^2 = -κ. This functional form provides a physical picture of the waiting time associated with (rare) transitions by factorizing the rate J into a kinetic part ∝ ωD(v^*) and the probability density P(v^*) of encountering a bubble with volume v^*, i.e., a configuration that relaxes to the vapor or the liquid phase with equal probability. In the following, we will compute the probability P(v^*) to find a bubble of critical size and derive an analytical expression for the diffusion constant D(v^*) needed in the CNT rate expression.

## Free Energy of Cavitation at Negative Pressures

Using umbrella sampling simulations, we have computed the equilibrium bubble density ρ(v) at a temperature T = 296.4 K and various pressures (Materials and Methods). For large bubbles, ρ(v) is equal to the probability density P(v)/V for the volume of the largest bubble as needed in Eq. 1 (42). The equilibrium bubble density is related to the Gibbs free energy g(v) of a bubble of volume v by $g(v) = -k_B T \ln[\rho(v)/\rho_0]$, where ρ_0 is a constant included to make the argument of the logarithm dimensionless. The value of ρ_0 is fixed by requiring that the Gibbs free energy of a bubble of size v = 0 vanishes. Note that the constant ρ_0, required to relate the cavitation free energy g(v) to the equilibrium bubble density ρ(v), is not specified in the framework of CNT. Various choices for ρ_0 have been made in the literature without

rigorous justification, as discussed in Materials and Methods. Here we use information from molecular simulations to determine the value of ρ_0 unambiguously (Supporting Information).

We obtain a quantitative description of the cavitation free energy within CNT by examining the free energetic cost of the bubble interface, i.e., the free energy without the mechanical work pv gained from expanding the system under tension, per surface area (Fig. S2). Remarkably, the free energetic cost of the vapor-liquid interface is independent of pressure within the accuracy of our computations (Fig. S3), and as such, for the wide range of pressures investigated, the free energy of cavitation differs only by the mechanical work pv. We find that CNT describes the free energy of bubble nucleation accurately, provided that the curvature dependence of the surface tension γ is taken into account. In particular, the free energy of cavitation is reproduced by

$$g(v) = 4\pi r^2(v) \frac{\gamma_0}{1 + 2\delta/r(v)} + pv, \quad [2]$$

where r(v) = (3v/4π)^1/3 is the radius of a sphere with volume v. Here, the parameters γ_0 = 20.24 k_B T/nm^2 and δ = 0.195 nm are obtained from a fit to the free energetic cost of the liquid-vapor interface. Bubble free energies g(v) for various pressures as well as the estimates from Eq. 2, which agree almost perfectly with the simulation data (dashed black lines), are shown in Fig. 1. Over the range of bubble volumes studied here, the value of δ obtained from the fit is positive, which indicates that the concave curvature of the interface decreases the surface tension γ, thereby favoring bubbles over droplets (a discussion of the curvature dependence of the surface tension is provided in Supporting Information).

## Bubble Morphology

At the conditions studied here, bubbles are essentially voids in the metastable liquid, which, for bubble volumes v ≲ 10 nm^3, rarely contain vapor molecules (34, 35). Visual inspection indicates that small bubbles mostly have an irregular shape, which becomes more compact as the bubbles grow larger (some representative bubbles of different size are depicted in Fig. 24). Larger bubbles are predominantly compact and may be viewed as resembling spheres with strongly undulating surfaces (34, 35). This observation is confirmed by computing the average asphericity of

<figure><img src="image_3.png" /><figcaption>Fig. 1. Free energy g(v) of bubble nucleation as a function of bubble volume v for various negative pressures at T = 296.4 K obtained from umbrella sampling calculations. Dashed lines indicate CNT predictions from Eq. 2, which describes the free energy very accurately over the investigated pressure range. In the framework of CNT, the curves can be understood as a result of the competition between the free energetic cost of forming the liquid-vapor interface (which dominates in the case of small bubbles) and the mechanical work gained from expanding the system under tension (favoring large bubbles). The location of the resulting maximum in the free energy corresponds to the volume of the critical bubble, v^*: bubbles of this volume are least likely to be encountered in an equilibrium configuration, and overcoming this free energy barrier is the rate-limiting step in cavitation away from the spinodal.</figcaption></figure>
bubbles defined as $\alpha = \lambda_{max}/\lambda_{min} - 1$, where $\lambda_{max}$ and $\lambda_{min}$ are the largest and smallest eigenvalues of the gyration tensor of the bubble, respectively. As shown in Fig. 2A, the asphericity is only weakly dependent on pressure and decreases with increasing bubble volume.

The free energetic cost of forming bubbles in water is inti- mately connected to breaking and rearranging hydrogen bonds (HBs) at the interface. The hydrogen bonding structure at the liquid-vapor interface depends on the size of the bubble (43, 44). For small bubbles, HBs in the liquid are rearranged, and the fraction of broken HBs at the interface is similar to that of the bulk liquid, whereas in the case of large bubbles, the bubble surface becomes similar to the flat vapor-liquid interface. As shown in Fig. 2B, the number of broken HBs per molecule at the interface increases with bubble size, and the fraction of free OH groups at the interface decays roughly linearly with its mean curvature $r^{-1}$ over the studied range of bubble volumes.

## Bubble Dynamics

Because CNT with a curvature-dependent surface tension describes the free energy of cavitation very accurately, thus providing the volume $v^*$ of the critical bubble and the curvature $-\omega^2$ of the barrier, all that is needed to predict rates via Eq. 1 is the diffusivity $D(v^*)$ of the bubble volume in the barrier region. In the following, we use the Rayleigh-Plesset (RP) equation (45-47), which describes the dynamics of a vapor bubble in a fluid at the macroscopic level, to derive an analytical expression that

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Shape and hydrogen bonding structure of bubbles. (A) Asphericity $\alpha$ as a function of bubble volume from configurations obtained via umbrella sampling. By construction, $\alpha$ is zero for a perfect sphere, and higher values indicate shapes with higher aspect ratios. Also shown are bubbles (not to scale) observed at $p = -150$ MPa whose asphericities and volumes are indicated by arrows. (B) Fraction of free OH groups $n_{OH}/n_{mol}$ at the bubble surface as a function of the inverse radius $r^{-1}$ of a sphere with volume $v$. The arrow indicates the fraction $n_{OH}/n_{mol}$ for a flat interface at 300 K at ambient pressure from ref. 64. Note that we give the fraction of broken hydrogen bonds per molecule, so the highest possible value is 2. The depicted configuration contains a bubble of critical size at a pressure of $p = -150$ MPa, where the yellow spheres indicate the unoccupied grid points forming the largest bubble.</figcaption></figure>

relates the microscopic diffusion constant $D(v^*)$ to the macroscopic properties of the liquid.

The RP equation is the equation of motion for the volume $v$ of a spherical bubble evolving with internal pressure $p_b$ in a liquid with mass density $m$, viscosity $\eta$, and surface tension $\gamma$:

$$m\ddot{v} - \frac{mv^2}{6v} = 4\pi \left(\frac{3v}{4\pi}\right)^{\frac{1}{3}} \left[ p_b - p - 2\gamma \left(\frac{4\pi}{3v}\right)^{\frac{1}{3}} - \frac{4\eta}{3} \frac{\dot{v}}{v} \right]. \quad [3]$$

Here, for simplicity, we neglect the curvature dependence of the surface tension but stress that the following derivation can be easily generalized (Supporting Information), and all results shown in Figs. 3 and 4 were obtained including this correction. Neglecting the inertial terms on the left-hand side of the RP equation, one finds

$$\dot{v} = -\frac{3v}{4\eta} \left[ p + 2\gamma \left(\frac{4\pi}{3v}\right)^{\frac{1}{3}} \right] = -\frac{1}{\Gamma(v)} \frac{dg(v)}{dv}, \quad [4]$$

where we assumed that the pressure inside the bubble is negligible. In the above equation we have rewritten the right-hand side to indicate that the time evolution of the volume $v$ can be viewed as an overdamped motion on the CNT free energy $g(v) = (36\pi v^2)^{1/3}\gamma + pv$ under the effect of the friction $\Gamma(v) = 4\eta/3v$.

Because thermal fluctuations play an important role for microscopic bubbles, the RP equation is augmented with a random force $F(t) = \sqrt{2k_B T/\Gamma(v)\xi(t)}$, where $\xi(t)$ is Gaussian white noise and the magnitude of the force is determined by the fluctuation-dissipation theorem. The diffusion coefficient for the bubble volume then follows from the Einstein relation, $D(v) = 3k_B T v/4\eta$ (note that this result also holds if the surface tension depends on the mean curvature of the bubble). Inserting the critical $v^* = 32\pi\gamma^3/3|p|^3$, we finally obtain the diffusion coefficient at the top of the barrier needed for the rate calculation

$$D(v^*) = \frac{8\pi k_B T \gamma^3}{\eta |p|^3}. \quad [5]$$

Including the curvature dependence of the surface tension for $v^*$ and $\gamma$ yields a similar but slightly more complicated formula (Supporting Information).

A comparison between the diffusion constant $D(v^*)$ obtained from the RP equation combined with CNT and the estimate obtained directly from simulation (Materials and Methods) is shown in Fig. 3. The viscosity at negative pressures needed in the formula for the diffusion constant was determined in molecular dynamics simulations using the Green-Kubo relation (Fig. S4). The analytical formula obtained from the RP-CNT approach underestimates the diffusivity only by about a factor of 2 compared with simulation results, which is remarkable considering that this estimate is obtained from a macroscopic approach based on hydrodynamics. Moreover, by virtue of the pressure dependence of $v^*$ in CNT, it predicts the scaling $\propto |p|^{-3}$ of the diffusion constant with pressure accurately, suggesting that the dynamics of bubble growth are essentially controlled by the viscosity of the liquid.

## Cavitation Rates

We are now in a position to predict cavitation rates according to Eq. 1 over a wide range of pressures, including the strongest tensions observed in experiment. As a point of comparison, we have computed cavitation rates numerically using a method akin to the divided-saddle method (48) based on the Bennett-Chandler (BC) (49, 50) approach and transition interface sampling (TIS) (51), respectively (Materials and Methods). The obtained cavitation rates, shown in Fig. 4, vary by more than 30 orders of
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** The diffusion constant $D(v^*)$ on top of the free energy barrier obtained from the Rayleigh-Plesset equation predicts the correct scaling with pressure $p$. The RP estimate (orange line) was obtained by using the volume $v^*$ of the critical bubble and the curvature $-w^2$ of the barrier from CNT, including a curvature dependent surface tension. The scaling behavior of the diffusion constant obtained from simulation (blue squares) is illustrated by a fit $\propto p^{-3}$ (dashed black line).</figcaption></figure>

magnitude over the studied range of pressures. The numerical results are accurately reproduced by CNT based on Eq. 1 with a curvature-dependent surface tension and the correct value of $\rho_0$ as well as the kinetic prefactor from the RP equation. In contrast, “plain” CNT, i.e., CNT with a constant surface tension and a commonly used expression for $\rho_0$ (*Materials and Methods*), underestimates the cavitation rates by more than 15 orders of magnitude. This shortcoming illustrates the importance of including microscopic information, such as a curvature-dependent surface tension and the correct value of $\rho_0$, for the accurate prediction of rates.

By computing the cavitation pressure $p_{cav}$ from the rates shown in Fig. 4 we can directly compare the results obtained here to the conflicting experimental estimates for the limit of metastability of water under tension. The obtained estimate for the cavitation pressure $p_{cav} \approx -126$ MPa is in line with the results obtained in inclusion experiments (1–6). In contrast, the predicted cavitation tension is more negative by about 100 MPa than the data obtained via other experimental techniques would suggest (20, 22). Because the simulation setup excludes impurities in the fluid by design, this suggests that cavitation in these cases is indeed heterogeneous as was suspected in previous works (4, 21), which explains the significantly lower stability of water under tension in these experiments (a detailed discussion is provided in *Supporting Information*).

## Conclusions

At ambient temperature and strong tension, bubbles in metastable water are essentially voids in the liquid whose shape can deviate significantly from the assumption of a spherical nucleus made in CNT, depending on their size. Nonetheless, provided the dependence of the surface tension on the average curvature is included, the free energetics of bubble formation can be quantitatively described in the framework of CNT. We find that the curvature contribution favors the cavity over the droplet, i.e., $\delta > 0$, in agreement with experimental results (4). In light of conflicting results on the sign of $\delta$ in water, further study is required to elucidate the influence of the chosen water model and biasing toward certain cavity shapes on the obtained value of $\delta$.

By including the effect of thermal fluctuations in the Rayleigh-Plesset equation, we obtain an estimate for the bubble diffusivity that accurately reproduces the pressure dependence found in simulation and scales inversely with the viscosity of the liquid. Combining the kinetic prefactor determined for this diffusivity with the equilibrium bubble density yields a CNT expression for the cavitation rate that reproduces the nucleation rates very well for negative pressures. However, the microscopic mechanism for cavitation is expected to change for higher pressures

and temperatures, where the saturated vapor density is significantly higher than at the temperature studied here. At those conditions, similarly to droplet nucleation (52), the transport of molecules across the interface via evaporation and condensation will have a stronger influence on the kinetics of bubble growth, thereby diminishing the influence of viscous damping on the dynamics of the bubble.

The estimate for the cavitation pressure obtained from our rate calculations agrees well with the data from inclusion experiments, thus calling the conflicting results harvested by other techniques into question. Because the latter methods greatly underestimate the stability of water under tension, heterogeneous cavitation due to impurities is a likely explanation for this discrepancy.

## Materials and Methods

**Simulation Details.** We simulate $N=2,000$ water molecules in the isothermal-isobaric ensemble at a temperature of $T=296.4$ K using a rigid, non-polarizable 4-site transferable intermolecular potential (TIP4P/2005) (53), where the long-range interactions are treated with Ewald summation. The rate computations are carried out using molecular dynamics by integrating the equations of motion with a time step of 2 fs using a time-reversible quaternion-based integrator that maintains the rigid geometry of water molecules (54). Constant pressure is ensured by a barostat based on the Andersen approach (55) coupled to a Nosé-Hoover thermostat chain (56). Equilibrium free energies are computed by use of umbrella sampling (US) in conjunction with the hybrid Monte Carlo (HMC) (57) scheme. Here we use a modified version of the Miller integrator (58) with a Liouville operator decomposition according to Omelyan (59), which reduces fluctuations in the total energy significantly, thereby allowing the use of a time step of 7 fs. Each HMC step consists of three MD integration steps, constant pressure was implemented by isotropic volume fluctuations according to the Metropolis criterion, and sampling was enhanced by replica exchange moves (60) between neighboring windows. For the direct computation of cavitation rates we use transition interface sampling (TIS) (51), where we implemented time reversal and replica exchange moves in addition to shooting moves (described in detail in refs. 41, 61, and 62). The probability histograms for the individual windows in US and TIS were spliced together using a self-consistent histogram method (63).

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Predictions obtained from CNT using microscopic information are in excellent agreement with cavitation rates $J$ from direct simulation. The estimates obtained from simulations by a variant of the Bennett-Chandler method (blue squares) agree well with the transition interface sampling (red circles) reference calculations (*Materials and Methods*). Predictions of curvature-corrected CNT (orange line) with the correct value of $\rho_0$ using the kinetic prefactor shown in Fig. 3 yield excellent agreement with simulation results, whereas plain CNT (gray line) severely underestimates the cavitation rate. For plain CNT, we chose $\rho_0 = n_l n_v$, where $n_l$ and $n_v$ are the number density of the liquid and the vapor, respectively (68). These rate estimates allow for a direct comparison with conflicting experimental predictions on the stability of water under tension by computing the cavitation pressure $p_{cav}$. Following ref. 22, we define $p_{cav}$ such that the probability to observe a cavitation event is $P=1/2$ in a system of volume $V = 1,000 \ \mu m^3$ over an observation time of $\tau = 1$ s. Assuming that the cavitation events are associated with an exponential waiting time, as is typical for activated processes, a rate of $J = \ln 2/(V\tau)$ (dashed black line) is compatible with this requirement. Its intersection with the CNT prediction gives the cavitation pressure $p_{cav} \approx -126$ MPa.</figcaption></figure>
**Order Parameter.** We study homogeneous bubble nucleation from overstretched metastable water using the volume of the largest bubble as a local order parameter. Estimates for the volume $v$ of each bubble present in the system are obtained by use of the V-method, which was developed to give thermodynamically consistent estimates for the bubble volume (35). (Note that the nomenclature was adapted to facilitate readability: $v/\xi$ in this work corresponds to $V/V_{\text{bubble}}/v$ in ref. 35.) The V-method is a grid-based clustering approach to bubble detection (30), calibrated such that its estimate $v$ for the volume of a bubble corresponds to the average change in system volume due to the presence of such a bubble:

$$v(\xi) = \frac{\partial}{\partial n} \langle V \rangle_{n(\xi)}. \qquad [6]$$

Here $\xi$ is the preliminary bubble volume estimate from the grid-based method, i.e., the total volume of all vapor-like grid cubes belonging to the bubble, and $\langle V \rangle_{n(\xi)}$ is the average volume of the system when $n$ bubbles of size $\xi$ are present. As such, $v(\xi)$ corresponds to the average change in system volume when a single bubble of size $\xi$ is added to or removed from the system. For large bubbles, i.e., for bubble volumes where $n(\xi)$ is either zero or one and there are no larger bubbles present in the system, Eq. 6 becomes

$$v(\xi) = \langle V \rangle_{\xi} - \langle V \rangle, \qquad [7]$$

where $\langle V \rangle_{\xi}$ is the average volume of the system when the largest bubble is of size $\xi$ and $\langle V \rangle$ is the average volume of the unconstrained metastable liquid at the thermodynamic state point.

On average, because the vapor density in the interior of bubbles is negligible, volume estimates obtained by Eq. 7 are equal to those obtained by computing the equimolar dividing surface between liquid and the largest cavity for each configuration. As a result, the obtained estimates for the bubble volume fulfill the nucleation theorem (23), i.e., $\partial g(v^*)/\partial p = v^*$, and $pv$ corresponds to the mechanical work gained with respect to the metastable liquid by expanding the system volume at negative pressures. Details on the calibration of the V-method for the state points investigated in this work are given in Fig. S5.

**Bubble Density.** To compute the equilibrium bubble density $\rho(v)$, we first carry out a straightforward molecular dynamics simulation and compute $\langle n(v, \Delta v) \rangle$, the average number of bubbles with a volume in a narrow interval $[v, v + \Delta v]$. To compute $n(v, \Delta v)$ for larger bubbles which do not form spontaneously on the time scale of the simulation, we carry out umbrella sampling simulations with a bias on the volume of the largest bubble. The resulting curves are joined, thus yielding $\rho(v) = \langle n(v, \Delta v) \rangle / \langle \langle V \rangle \Delta v \rangle$ over a wide range of bubble volumes.

**Detecting Hydrogen Bonds at the Liquid-Vapor Interface.** We identify molecules as belonging to the bubble surface when they are within 3.5 Å of the bubble. This cutoff radius is identical to the radius of the exclusion spheres used to determine occupied grid points during the evaluation of the order parameter (for an in-depth description, see ref. 35), and thus, all water molecules forming the boundary layer in our bubble detection procedure are part of the interface. When analyzing whether two water molecules form a hydrogen bond with each other, we use the criterion used in ref. 64 in a study of the flat vapor-liquid interface to facilitate easy comparison between the obtained results. For molecule A to be considered as donating a hydrogen bond to molecule B, two criteria have to be fulfilled simultaneously: The distance between the oxygens $d_{O_A O_B} < 3.5$ Å and the maximum angle $\text{O}_A - \text{H} \cdots \text{O}_B > 140^\circ$.

**Rate Calculation.** We use a method based on the Bennett-Chandler approach (49, 50) to obtain rates estimates without any assumptions about the dynamics of the bubble in the liquid. In addition to the states A (metastable liquid) and B (far enough to the right of the free energy barrier such that the system is committed to transitioning to the vapor phase), we introduce a state S around the dividing surface, akin to the approach taken in the divided-saddle method (48). An ensemble of trajectories, each L steps long, is generated by propagating checkpoints selected from the region S forward and backward in time. From these trajectories one then computes the time

correlation function $C_{AB}(t)$, which is the conditional probability to find the system in B at time t provided it is in A at time 0,

$$C_{AB}(t) = (L+1) \frac{\langle h_A(0)h_B(t) \rangle}{N_S[\langle x \rangle]} \langle h_S \rangle_G \langle h_A \rangle_G. \qquad [8]$$

Here $h_{A/B}$ is 1 when the system is in state A/B and zero else, $N_S[\langle x \rangle]$ is the number of configurations of a trajectory $x(\tau)$ in the saddle domain, and $\langle \cdots \rangle_G$ denotes an average over the trajectories generated from points in S. The ratio $\langle h_S \rangle_G / \langle h_A \rangle_G$ is the equilibrium probability of finding the system in S relative to the equilibrium probability of state A, and it can be determined from the free energy $g(v)$. The transition rate constant $k_{AB}$ is then obtained by computing the numerical derivative $dC_{AB}/dt$ in the time range where $C_{AB}(t)$ is linear.

Nucleation rates calculated at $p = -165$ MPa and $-150$ MPa using transition interface sampling (51) (TIS; red circles in Fig. 4) agree with the estimates of the BC-based approach up to statistical error. As an additional point of comparison, we used the BC-based approach to compute rates at $T = 280$ K and $p = -225$ MPa, where nucleation is spontaneous on the time scale of an unconstrained molecular dynamics simulation starting in the metastable liquid. The estimate $J = 3.1 \times 10^{-5} \text{ps}^{-1} \cdot \text{nm}^{-3}$ obtained from straightforward MD simulations in ref. 34 agrees well with the BC-based estimate of $J = 7.4 \times 10^{-5} \text{ps}^{-1} \cdot \text{nm}^{-3}$.

**Computation of the Diffusion Constant.** Because the volume of the largest bubble is a good reaction coordinate for the transition, its diffusivity can be computed via mean first passage times (65, 66). Assuming that the diffusion coefficient does not change significantly in the barrier region, i.e., $D(v) = D(v^*)$, to second order it can be expressed as $D = b^2(1 - \beta b^2 \omega^2 / 6) / (2 \langle \tau(b) \rangle)$, where b is the distance of the absorbing boundary from the top of the free energy barrier, approximated by an inverted parabola with curvature $-\omega^2$, and $\langle \tau(b) \rangle$ is the mean first passage time for a given value of b. As a starting point at the top of the barrier we used equilibrium configurations created by umbrella sampling where the system contained a cluster of critical size and drew the particle velocities as well as the thermostat and barostat velocities at random from the appropriate Maxwell-Boltzmann distributions.

**Plain CNT.** As a point of comparison, we obtain an estimate for the cavitation rates from CNT with a constant surface tension $\gamma_0 = 17.09 \text{ k}_B T/\text{nm}^2$ for TIP4P/2005 water (67). The CNT estimate for the rate is given by

$$J = \frac{\sqrt{k_B T \gamma_0^3}}{\eta |p|} \rho_0 e^{-\beta 16 \pi \gamma_0^3 / 3p^2}. \qquad [9]$$

The equation above was obtained from Eqs. 1 and 5, where $\omega = p^2 / \sqrt{32\pi\gamma_0^3}$ and the probability density $P(v)/V = \rho_0 \exp(-\beta g(v^*))$. Here $g(v^*) = 16\pi\gamma_0^3 / 3p^2$, and the normalization constant was chosen as $\rho_0 = n_l n_v \approx 4.4 \times 10^{-3} \text{nm}^{-6}$, where $n_l$ and $n_v$ are the number density of the metastable liquid and the number density of the vapor at coexistence (68), respectively. Note that the prefactor $\rho_0$ is not uniquely defined in the framework of CNT, and various choices have been used in the literature (25, 26, 68). These choices lead to estimates ranging from $\rho_0 = 9.4 \times 10^{-14} \text{nm}^{-6}$ to $\rho_0 = 2.4 \times 10^8 \text{nm}^{-6}$ at $p = -135$ MPa (we obtain $\rho_0 = 0.02 \text{ nm}^{-6}$ from the simulation data shown in Fig. S2). The resulting predictions for the cavitation rates underestimate the values determined from simulation by 6–27 orders of magnitude.

**ACKNOWLEDGMENTS.** We thank S. Garde, P. Geissler, V. Molinero, A. Patel, E. Sanz, A. Tröster, E. Vanden-Eijnden, C. Vega, and S. Venkatari for insightful comments. Calculations were carried out on the Vienna Scientific Cluster. The work of G.M., P.G., and C.D. was supported by the Austrian Science Foundation (FWF) under Grant P24681-N20 and within the Spezialforschungsbereich Vienna Computational Materials Laboratory (Grant F41). P.G. also acknowledges financial support from FWF Grant P22087-N16, and F.C. acknowledges support from the European Research Council under the European Union's Seventh Framework Programme for Research and Technological Development Grant Agreement 240113. C.V. acknowledges financial support from a Marie Curie Integration Grant 322326-COSAAC-FP7-PEOPLE-CIG-2012 and a Ramon y Cajal tenure track. The team at Madrid acknowledges funding from the Ministerio de Economía y Competitividad Spanish Grants FIS2013-43209-P, FIS2016-78117-P, and FIS2016-78847-P.

1. Green JL, Durben DJ, Wolf GH, Angell CA (1990) Water and solutions at negative pressure: Raman spectroscopic study to -80 megapascals. *Science* 249(4969):649–652.

2. Zheng Q, Durben DJ, Wolf GH, Angell CA (1991) Liquids at large negative pressures: water at the homogeneous nucleation limit. *Science* 254(5033):829–832.

3. Alvarenga AD, Grimsditch M, Bodnar RJ (1993) Elastic properties of water under negative pressures. *J Chem Phys* 98(11):8392–8396.

4. Azouzi MEM, Ramboz C, Lenain JF, Caupin F (2013) A coherent picture of water at extreme negative pressure. *Nat Phys* 9(1):38–41.

5. Mercury L, Shmuloovich K (2014) Experimental superheating and cavitation of water and solutions at spinodal-like negative pressures. *Transport and Reactivity of Solutions in Confinement*, NATO Science for Peace and Security Series C: Environmental Security, eds Mercury L, Tas N, Zilberbrand M (Springer, Dordrecht, The Netherlands), pp 159–171.
6. Pallares G, et al. (2014) Anomalies in bulk supercooled water at negative pressure. *Proc Natl Acad Sci USA* 111(22):7936-7941.

7. Stroock AD, Pagay VV, Zwieniecki MA, Michele Holbrook N (2014) The physicochemical hydrodynamics of vascular plants. *Annu Rev Fluid Mech* 46(1):615-642.

8. Ponomarenko A, et al. (2014) Ultrasonic emissions reveal individual cavitation bubbles in water-stressed wood. *J R Soc Interface* 11(99):2014080.

9. Larter M, et al. (2015) Extreme aridity pushes trees to their physical limits. *Plant Physiol* 168(3):804-807.

10. Rowland L, et al. (2015) Death from drought in tropical forests is triggered by hydraulics not carbon starvation. *Nature* 528(7580):119-122.

11. Wheeler TD, Stroock AD (2008) The transpiration of water at negative pressures in a synthetic tree. *Nature* 455(7210):208-212.

12. Vincent O, Marmottant P, Quinto-Su PA, Ohl CD (2012) Birth and growth of cavitation bubbles within water under tension confined in a simple synthetic tree. *Phys Rev Lett* 108(18):184502.

13. Noblin X, et al. (2012) The fern sporangium: A unique catapult. *Science* 335(6074):1322.

14. Ohl CD, et al. (2006) Sonoporation from jetting cavitation bubbles. *Biophys J* 91(11):4285-4295.

15. Adhikari U, Goliaei A, Berkowitz ML (2015) Mechanism of membrane poration by shock wave induced nanobubble collapse: A molecular dynamics study. *J Phys Chem B* 119(20):6225-6234.

16. Ohsaka K, Trinh EH (1998) Dynamic nucleation of ice induced by a single stable cavitation bubble. *Appl Phys Lett* 73(1):129-131.

17. Yu D, Liu B, Wang B (2012) The effect of ultrasonic waves on the nucleation of pure water and degassed water. *Ultrason Sonochem* 19(3):459-463.

18. Kumar P, Saini R (2010) Study of cavitation in hydro turbines—A review. *Renew Sustain Energy Rev* 14(1):374-383.

19. Debenedetti PG (2013) Physics of water stretched to the limit. *Nat Phys* 9(1):7-8.

20. Herbert E, Balibar S, Caupin F (2006) Cavitation pressure in water. *Phys Rev E Stat Nonlin Soft Matter Phys* 74(4 Pt 1):041603.

21. Davitt K, Arvengas A, Caupin F (2010) Water at the cavitation limit: Density of the metastable liquid and size of the critical bubble. *EPL* 90(1):16002.

22. Caupin F, Herbert E (2006) Cavitation in water: A review. *C R Phys* 7(9-10):1000-1017.

23. Kashchiev D (1982) On the relation between nucleation work, nucleus size, and nucleation rate. *J Chem Phys* 76(10):5098-5102.

24. Caupin F (2005) Liquid-vapor interface, cavitation, and the phase diagram of water. *Phys Rev E Stat Nonlin Soft Matter Phys* 71(5 Pt 1):051605.

25. Zeng XC, Oxtoby DW (1991) Gas-liquid nucleation in Lennard-Jones fluids. *J Chem Phys* 94(6):4472-4478.

26. Oxtoby DW (1992) Homogeneous nucleation: Theory and experiment. *J Phys Condens Matter* 4(38):7627-7650.

27. Shen VK, Debenedetti PG (1999) A computational study of homogeneous liquid-vapor nucleation in the Lennard-Jones fluid. *J Chem Phys* 111(8):3581-3589.

28. Vishnyakov A, Debenedetti PG, Neimark AV (2000) Statistical geometry of cavities in a metastable confined fluid. *Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics* 62(1 Pt A):538-544.

29. Neimark AV, Vishnyakov A (2005) The birth of a bubble: A molecular simulation study. *J Chem Phys* 122(5):54707.

30. Wang ZJ, Valeriani C, Frenkel D (2009) Homogeneous bubble nucleation driven by local hot spots: A molecular dynamics study. *J Phys Chem B* 113(12):3776-3784.

31. Baidakov VG, Bobrov KS, Teterin AS (2011) Cavitation and crystallization in a metastable Lennard-Jones liquid at negative pressures and low temperatures. *J Chem Phys* 135(5):054512.

32. Meadly SL, Escobedo FA (2012) Thermodynamics and kinetics of bubble nucleation: Simulation methodology. *J Chem Phys* 137(7):074109.

33. Torabi K, Corti DS (2013) Toward a molecular theory of homogeneous bubble nucleation: II. Calculation of the number density of critical nuclei and the rate of nucleation. *J Phys Chem B* 117(41):12491-12504.

34. Abascal JLF, Gonzalez MA, Aragones JL, Valeriani C (2013) Homogeneous bubble nucleation in water at negative pressure: A Voronoi polyhedra analysis. *J Chem Phys* 138(8):084508.

35. Gonzalez MA, et al. (2014) Detecting vapour bubbles in simulations of metastable water. *J Chem Phys* 141(18):18C511.

36. Gonzalez MA, Abascal JLF, Valeriani C, Bresme F (2015) Bubble nucleation in simple and molecular liquids via the largest spherical cavity method. *J Chem Phys* 142(15):154903.

37. Zahn D (2004) How does water boil? *Phys Rev Lett* 93(22):227801.

38. Cho WJ, et al. (2014) Limit of metastability for liquid and vapor phases of water. *Phys Rev Lett* 112(15):157802.

39. Kramers H (1940) Brownian motion in a field of force and the diffusion model of chemical reactions. *Physica* 7(4):284-304.

40. Schulten K, Schulten Z, Szabo A (1981) Dynamics of reactions involving diffusive barrier crossing. *J Chem Phys* 74(8):4426-4432.

41. Dellago C, Bolhuis PG, Geissler PL (2002) Transition path sampling. *Advances in Chemical Physics* (John Wiley, New York), Vol 123, pp 1-78.

42. Maibaum L (2008) Comment on "Elucidating the mechanism of nucleation near the gas-liquid spinodal". *Phys Rev Lett* 101(1):019601, discussion 019602.

43. Lum K, Chandler D, Weeks JD (1999) Hydrophobicity at small and large length scales. *J Phys Chem B* 103(22):4570-4577.

44. Chandler D (2005) Interfaces and the driving force of hydrophobic assembly. *Nature* 437(7059):640-647.

45. Plesset MS, Prosperetti A (1977) Bubble dynamics and cavitation. *Annu Rev Fluid Mech* 9:145-185.

46. Kagan Y (1960) The kinetics of boiling of a pure liquid. *Russ J Phys Chem* 34(1):42-48.

47. Leighton TG (2008) The Rayleigh-Plesset equation in terms of volume with explicit shear losses. *Ultrasonics* 48(2):85-90.

48. Daru J, Stirling A (2014) Divided saddle theory: A new idea for rate constant calculation. *J Chem Theory Comput* 10(3):1121-1127.

49. Bennett CH (1977) Molecular dynamics and transition state theory: The simulation of infrequent events. *Algorithms for Chemical Computations* (Am Chem Soc, Washington, DC), pp 63-97.

50. Chandler D (1978) Statistical mechanics of isomerization dynamics in liquids and the transition state approximation. *J Chem Phys* 68(6):2959-2970.

51. van Erp TS, Moroni D, Bolhuis PG (2003) A novel path sampling method for the calculation of rate constants. *J Chem Phys* 118(17):7762-7774.

52. Becker R, Döring W (1935) Kinetic treatment of grain-formation in super-saturated vapours. *Ann Phys* 416(8):719-752.

53. Abascal JLF, Vega C (2005) A general purpose model for the condensed phases of water: TIP4P/2005. *J Chem Phys* 123(23):234505.

54. Kamberaj H, Low RJ, Neal MP (2005) Time reversible and symplectic integrators for molecular dynamics simulations of rigid molecules. *J Chem Phys* 122(22):224114.

55. Andersen HC (1980) Molecular dynamics simulations at constant pressure and/or temperature. *J Chem Phys* 72(4):2384-2393.

56. Tuckerman ME (2010) *Statistical Mechanics: Theory and Molecular Simulation* (Oxford University Press, Oxford).

57. Duane S, Kennedy A, Pendleton B, Roweth D (1987) Hybrid Monte Carlo. *Phys Rev Lett* 195(2):216-222.

58. Miller T, III, et al. (2002) Symplectic quaternion scheme for biophysical molecular dynamics. *J Chem Phys* 116(20):8649-8659.

59. Omelyan IP, Mryglod IM, Folk R (2002) Optimized Verlet-like algorithms for molecular dynamics simulations. *Phys Rev E Stat Nonlin Soft Matter Phys* 65(5 Pt 2):056706.

60. Geyer CJ, Thompson EA (1995) Annealing markov chain monte carlo with applications to ancestral inference. *J Am Stat Assoc* 90(431):909-920.

61. van Erp TS (2007) Reaction rate calculation by parallel path swapping. *Phys Rev Lett* 98(26):268301.

62. Bolhuis PG (2008) Rare events via multiple reaction channels sampled by path replica exchange. *J Chem Phys* 129(11):114108.

63. Ferrenberg AM, Swendsen RH (1989) Optimized Monte Carlo data analysis. *Phys Rev Lett* 63(12):1195-1198.

64. Vila Verde A, Bolhuis PG, Campen RK (2012) Statics and dynamics of free and hydrogen-bonded OH groups at the air/water interface. *J Phys Chem B* 116(31):9467-9481.

65. Berezhkovskii AM, Szabo A (2013) Diffusion along the splitting/commitment probability reaction coordinate. *J Phys Chem B* 117(42):13115-13119.

66. Lu J, Vanden-Eijnden E (2014) Exact dynamical coarse-graining without time-scale separation. *J Chem Phys* 141(4):044109.

67. Vega C, de Miguel E (2007) Surface tension of the most popular models of water by using the test-area simulation method. *J Chem Phys* 126(15):154707.

68. Blander M, Katz JL (1975) Bubble nucleation in liquids. *AIChE J* 21(5):833-848.

69. Tolman RC (1949) The effect of droplet size on surface tension. *J Chem Phys* 17(3):333-337.

70. Rowlinson JS, Widom B (1989) *Molecular Theory of Capillarity* (Dover, New York).

71. Tröster A, Oettel M, Block B, Virnau P, Binder K (2012) Numerical approaches to determine the interface tension of curved interfaces from free energy calculations. *J Chem Phys* 136(6):064709.

72. Joswiak MN, Duff N, Doherty MF, Peters B (2013) Size-dependent surface free energy and Tolman-corrected droplet nucleation of tip4p/2005 water. *J Phys Chem Lett* 4(24):4267-4272.

73. Bruot N, Caupin F (2016) Curvature dependence of the liquid-vapor surface tension beyond the tolman approximation. *Phys Rev Lett* 116(5):056102.

74. Sedlmeier F, Netz RR (2012) The spontaneous curvature of the water-hydrophobe interface. *J Chem Phys* 137(13):135102.

75. Vaikuntanathan S, Geissler PL (2014) Putting water on a lattice: The importance of long wavelength density fluctuations in theories of hydrophobic and interfacial phenomena. *Phys Rev Lett* 112(2):020603.

76. Factorovich MH, Molinero V, Scherlis DA (2014) Vapor pressure of water nanodroplets. *J Am Chem Soc* 136(12):4508-4514.

77. Wilhelmsen Ø, Bedeaux D, Reguera D (2015) Communication: Tolman length and rigidity constants of water and their role in nucleation. *J Chem Phys* 142(17):171103.

78. Lau GV, Hunt PA, Müller EA, Jackson G, Ford IJ (2015) Water droplet excess free energy determined by cluster mitosis using guided molecular dynamics. *J Chem Phys* 143(24):244709.

79. Nevins D, Spera FJ (2007) Accurate computation of shear viscosity from equilibrium molecular dynamics simulations. *Mol Simul* 33(15):1261-1266.

80. González MA, Abascal JLF (2010) The shear viscosity of rigid water models. *J Chem Phys* 132(9):096101.

81. Caupin F (2015) Escaping the no man's land: Recent experiments on metastable liquid water. *J Non-Cryst Solids* 407:441-448.

82. Pallares G, Gonzalez MA, Abascal JLF, Valeriani C, Caupin F (2016) Equation of state for water and its line of density maxima down to -120 MPa. *Phys Chem Chem Phys* 18(8):5896-5900.

83. Stan CA, et al. (2016) Negative pressures and spallation in water drops subjected to nanosecond shock waves. *J Phys Chem Lett* 7(11):2055-2062.