Water cavitation from ambient to high temperatures

Predicting cavitation has proved a formidable task, particularly for water. Despite the experimental difficulty of controlling the sample purity, there is nowadays substantial consensus on the remarkable tensile strength of water, on the order of −120 MPa at ambient conditions. Recent progress significantly advanced our predictive capability which, however, still considerably depends on elaborate fitting procedures based on the input of external data. Here a self-contained model is discussed which is shown able to accurately reproduce cavitation data for water over the most extended range of temperatures for which accurate experiments are available. The computations are based on a diffuse interface model which, as only inputs, requires a reliable equation of state for the bulk free energy and the interfacial tension. A rare event technique, namely the string method, is used to evaluate the free-energy barrier as the base for determining the nucleation rate and the cavitation pressure. The data allow discussing the role of the Tolman length in determining the nucleation barrier, confirming that, when the size of the cavitation nuclei exceed the thickness of the interfacial layer, the Tolman correction effectively improves the predictions of the plain Classical Nucleation Theory.


Scientific Reports
| (2021) 11:20801 | https://doi.org/10.1038/s41598-021-99863-z www.nature.com/scientificreports/ conditions, where the bubble radius is sufficiently larger than the thickness of the interfacial layer, CNT provides qualitatively correct results. However, the CNT nucleation barrier is found to be overestimated leading to large discrepancies with the experimentally observed cavitation rates and cavitation pressures 6,10 .
Recently, Menzl et al. 11 corrected the CNT by including the Tolman length δ in the basic model to account for the effect of interface curvature on the surface tension σ . Assuming δ > 0 , as generally accepted, the curvature decreases σ and the nucleation free-energy barrier, �� * = 16/3πσ 3 * /�p 2 , is decreased. Since the barrier height controls the escape rate from the potential well, an increase of the nucleation rate by orders of magnitude follows.
In 11 , the classical Kramers approach was complemented by determining the diffusion coefficient in the Langevin equation for the bubble volume on the basis of the overdamped Rayleigh-Plesset equation for the nucleating bubble. This extended approach, completed with data from molecular dynamics (MD) simulations proved successful in reproducing the order of magnitude of the cavitation pressure estimated from quartz inclusion experiments at ambient temperature 5 .
In the present paper we attack the problem of homogeneous vapor bubble nucleation in water with a Diffuse Interface (DI) approach 12 that we recently developed for model fluids (Van der Walls equation of state or Lennard-Jones fluids) [13][14][15][16][17] . Historically, the DI model originated from the pioneering work on capillarity of Van der Waals who introduced a continuous, sharply varying density distribution ρ(r) in place of a sharp interface. This approach is a simplified version of density functional theory (DFT). The purpose of the proposed model is to reproduce the two crucial observables related to water cavitation, namely the nucleation rate J, defined as the number of bubbles nucleated per unit time and volume, and the cavitation pressure p cav . Given the volume of the liquid sample, V l , the latter is defined as the liquid pressure for which the probability to observe at least one cavitation bubble in a given time window, 0 ≥ t ≥ τ , equals 1/2. This definition is based on assuming nucleation as a random Poisson point process 18 , where single nucleation events are independently distributed in space and time.
The DI model is completed with a realistic equation of state for water-the IAPWS-95 EoS 19 -and the surface energy of the planar interface, σ 0 (T) , as a function of temperature 20 . The model implicitly includes the effect of interface curvature on the surface energy, accounting for the Tolman length, δ (see "Methods").

Results
One of the main results of the paper is provided in Fig. 1 which shows the cavitation pressure of ultra-pure water as a function of temperature. Introducing the nucleation rate J, the average number of nucleated bubbles observed in a liquid volume V l within the time window 0−−τ is, by definition, �n� = JV l τ . Given the Poisson distribution of the bubble number π(n|V l , τ ) = 1/n!(JV l τ ) n exp(−JV l τ ) , the probability of having at least one bubble is 1 − exp(−JV l τ ) . Based on the above assumptions, the cavitation pressure is the liquid pressure at which this probability is exactly 1/2, i.e. J(p cav ) = ln(2)/(V l τ ).
The data shown by the red thick curve in the figure were obtained by evaluating, for given temperature, the nucleation rate-see "Methods"-for different pressures to extract p cav by interpolation. Within the shaded region, obtained by varying V l τ in the range 1 ÷ 10 6 µm 3 · s , the cavitation pressures changes by less than 10% for all the analyzed temperatures, confirming the robustness of this observable.
The predictions of the present model compare favorably with the experimental data at high temperature 2 (quartz inclusions) 21 , (heat pulse method), see also the figure inset. As the temperature is reduced, the experimental scatter on the cavitation pressure data increases considerably. The point indicated by the blue symbol on Figure 1. Cavitation pressure as a function of temperature. The red solid line reports the values as predicted by the diffuse interface model, as explained in the text; the light red band shows the sensitivity to the experimental parameter V l τ = 1000 µm 3 · s when changed by a factor 10 −3 (lower limit) and 10 3 (upper limit). The dashdotted black, and the dashed blue, lines show the spinodal and saturation pressures, respectively, as provided by IAPWS EoS. The different symbols (green stars 2 , black dots 21 ) represent available experimental data. The arrow indicates the correction due to the matrix compliance effect. The two dotted curves report data from 22 exploiting two different water EoS, as indicated in the legend. The inset zooms into the high temperature region, close to the critical point.  Figure 2 provides the liquid mass density at cavitation as a function of temperature, namely the density of the liquid far away from the bubble (i.e. in the bulk liquid state corresponding to the assigned temperature and pressure). The numerical data shown by the red line correspond to those given in terms of pressure in the previous Fig. 1. The experimental data come from four different experiments. The plus symbols 23 provide cavitation density data from quartz inclusion experiments. They substantially correspond to the pressure data denoted by green stars 2 in the pressure diagram. The left-most point is presumably affected by some inaccuracy 2 and it is expected to underestimate the pressure, see the vertical green arrow in Fig. 1. The remaining experimental data points, see caption, provide data from acoustic experiments, which are believed to anticipate the cavitation. The inset shows the vapor density at the bubble center, an information which is unaccessible to experiments. Our model shows that vapor is close to saturation conditions, confirming expectations and CNT assumptions.
As explained in "Methods", the string method is used to find the critical state. Briefly, the string method 24,25 is a specialized approach to extract the minimum free energy path in the transition between (meta)stable states and allows to determine the free-energy barrier. This class of algorithms is suggested by the rare event nature of the transition process which, under certain circumstances, can be extremely slow, hindering the adoption of direct simulation strategies. The nucleation free-energy barriers for different temperatures in the range 50-350 • C are plotted in Fig. 3 as a function of the metastability level, µ lev = (µ − µ sat )/(µ sp − µ sat ) 26 , with µ the chemical potential and the subscripts sp and sat denoting spinodal and saturation conditions, respectively. The * vs µ lev curves tend to collapse, although not exactly, at high temperature, a feature consistent with the density functional theory (DFT) results reported in 26 for a Lennard-Jones fluid. In the present case, this behavior is apparently violated when the temperature decreases, red curve. The symbols in the main plot identify the barrier at the cavitation pressure, i.e. µ lev = (µ(p cav , T) − µ sat (T))/(µ sp (T) − µ sat (T)) . Clearly, the barrier in thermal units at cavitation conditions is substantially independent of temperature, owing to the definition of cavitation pressure ( J(p cav ) = ln(2)/(V l τ ) ) and the expression of the nucleation rate J = Ŵ 0 exp(−�� * /(k B T)) , which implies that the temperature dependence of * is given by the additive contribution ln(Ŵ 0 ) , where Ŵ 0 (T) is, e.g., the Blander and Katz 27 prefactor.
In a famous paper, Kashchiev 28 derived the celebrated nucleation theorem, (1/m)∂�� * /∂�µ = −n * , where n * is the molecule number in the critical cluster, with m the molecule mass, stating that it is expected to hold for classical, atomistic, homo-or heterogeneous, three-or two-dimensional nucleation. A more general result has been obtained in 29 where the theorem has been proved in the context of DFT, with continuously varying density fields. In 26 the original expression was modified by substituting the molecule number with the excess/defect of molecules, n * , to extend its range of application to bubbles.
The nucleation theorem is shown to be preserved also in the present DI context, as shown in the inset of  For each temperature, the data shown in the Fig. 4, symbols, are taken from the critical bubble at different metastability level. The lines are the fitting of the surface tension using the Tolman expression, σ (T, R) = σ 0 /(1 + 2δ/R) , where the parameter δ is the Tolman length 31 and σ 0 is the surface tension of the planar interface. We stress that σ 0 at the different temperatures are taken from the fitted IAPWS expression for water 20 . This uniquely determines the capillary parameter entering the DI free energy, (T) , Methods. Apparently,  Beside the cavitation pressure, the nucleation rate J is the other most important observable. The classical Kramers' theory 32 provides the reaction rate by modeling the process as the escape of a random walker from the free-energy potential well. Figure 5 shows J as obtained from the barrier heights taken from the present DI model using the prefactor Ŵ 0 = (ρ L ρ V /m 2 ) k B Tσ 3 /(η�p) . The solid lines in Fig. 5 show, for comparison, the data obtained using the Tolman-corrected CNT 11 while the broken lines correspond to the plain CNT.

Discussion
The diffuse interface method adopted in the present paper assumes the free energy of the capillary system as described by a functional of the density field ρ(r) expressed as a volume integral of the bulk free-energy density and a capillary contribution proportional to the squared density gradient, see "Methods". The equilibrium density field presents a strong gradient, sharply raising from the vapor to the liquid density at the interface between the two phases. The model can be interpreted as a gradient expansion of the functional used in DFT 33 . In can be considered as an effective theory where a few basic ingredients need to be included, namely the equilibrium thermodynamic properties of the liquid/vapor system. This allows the flexibility required to model realistic fluids, in the present case water, as described by the IAPWS EoS 19 for bulk free energy and surface tension, in a wide range of temperatures and pressures. In the model the EoS should encompass all the states along the transition between the two phases (from the metastable liquid to the stable vapor in the cavitation bubble). In particular, also the unstable states below the spinodal line need to be considered. This raises an issue, since no direct experimental information is available in the unstable region. As discussed in the "Methods", the procedure we have adopted consists in extrapolating data from the EoS to spinodal conditions and connect the limiting states across the unstable region. The interpolation is achieved with the constraints that (i) the pressure should be a strictly decreasing function of density at constant temperature (to prevent unphysical regions of local stability below the spinodal line), (ii) ∂p/∂ρ(ρ sp )| T = 0 (as required for spinodal states), and (iii) enforcing the continuity of pressure and free-energy. This procedure leaves a large freedom on the choice of the specific interpolation function. As shown in the Supplementary Information (SI), the model is insensitive to the details in the unstable region. The sensitiveness is larger to modifications of the shape of the spinodal curves, SI Fig. 2. In 22 , using an approach similar to the present one, the effect on cavitation of the topology of the spinodal line was discussed in some detail, concluding that a so-called reentrant spinodal, as in the case of the IAPWS EoS, is presumably better suited to capture cavitation data.
As shown in the "Results" section, our simulations reproduce remarkably well the best data available for bulk water. In comparing simulation data with experiments, one should be aware of certain technical difficulties encountered on the experimental side. As the authors of the experiments themselves underline, in certain cases, particularly at low temperature, it turned out difficult to strictly control the purity of the water sample. This is especially true of experiments done with acoustic excitations 3,4 . Water inclusions in quartz 1,2 , on the other hand, reduce the probability that impurities are present in the sample by using an extremely small volume of water, exploiting the quartz hydrophilicity to prevent bubble nucleation at the solid interface, which would lower the transition barrier. Nevertheless, each particular sample showed a different behavior in terms of cavitation pressure. Substantial effort was put in identifying the most suitable sample to perform repeated nucleation experiments in order to extract statistically reliable data 5 . Based on these considerations, the experimental data should, in general, be interpreted as overestimating the actual cavitation pressure pointing the attention to the low Figure 5. Bubble nucleation rate J as a function of the metastable liquid pressure for different temperatures. The DI data are reported with symbols, while the dashed and solid curves correspond to the predictions of the plain-CNT and the Tolman-corrected-CNT, respectively. The horizontal black line indicates the nucleation target rate J 1/2 = ln2/(V l τ ) , corresponding to the value of the rate such that one bubble is observed in a system of volume V l over an observation time τ with probability 1/2, following Ref. 6 . The black light band shows the sensitivity to the experimental parameters V l τ = 1000µm 3 × 1s by a factor of 0.001 (upper band) and 1000 (lower band). www.nature.com/scientificreports/ pressure envelope of the available data. With these comments in mind, the present simulation data are entirely consistent with experiments throughout the extended range of temperatures from 25 to 350 • C where reliable experimental data are available. As discussed in the introductory material, Kramers' theory was revised and adapted to bubble nucleation in a recent paper 11 where the Tolman length and the flat interface surface tension were fitted on MD data obtained with the TIP4P water model, predicting the cavitation pressure of −126 MPa at the temperature T = 23.25 • C . We like to stress that the Tolman length δ is implicitly included in the present DI model, where the density distribution ρ(r) is influenced by the size, hence the curvature, of the vapor bubble. As a consequence the excess grand potential that controls the actual surface tension felt by the system depends on bubble curvature. Figure 4 shows that the surface tension of the critical nucleus can be fitted remarkably well by the Tolman expression (see the caption for the values of δ ), at least in conditions where the interfacial thickness remains sufficiently small with respect to the bubble size for the very concept of Tolman length to make sense. As a consequence, at the temperature of 25 • C our model predicts p cav = −127.6 MPa without introducing additional information apart from the EoSs, Fig. 1.
The thickness ℓ of the interfacial layer is amenable to direct measurement with X-ray scattering and ellipsometry 34,35 respectively. The related data obtained from the present model are listed in the SI, Table 2, where the possibility to alter the interfacial thickness for given surface tension by changing the EoS in the unstable region is also briefly addressed. In fact, the interfacial thickness may easily be changed by 15% with no alteration of the nucleation rate and cavitation pressure.
As a final remark, we stress that several different prefactors, see, e.g. 11,27 , may be considered in the expression of the cavitation rate. The cavitation pressure is confirmed to be almost independent of the specific prefactor (SI, Fig. 3).

Conclusions and perspectives
Modeling cavitation in water has longly been an elusive issue. In this paper we have shown in detail how a diffuse interface method completed with realistic equations of state for bulk water and for the vapor-liquid surface tension can be exploited in combination with rare event techniques to evaluate the free-energy barrier.
The model is shown able to capture the experimental values of the cavitation pressure in water over the wider range of temperatures over which experimental data are available, e.g. from 25 • C to critical conditions. Our data provide a lower envelope for the low temperature data, which are known to overestimate the cavitation pressure and the particularly accurate low temperature (ambient) data by 5 are well matched by the present predictions.
We also confirm the findings discussed in 11 where the authors show that the effect of interfacial curvature is crucial to reproduce realistic nucleation barriers. We stress, however, that the Tolman parameter is not a fitting degree of freedom in the present model, which is completely determined by equation of state and planar interface free-energy. The present data allow to directly verify the validity of the Tolman correction which, as shown in Fig. (4), works perfectly well for larger sizes, deteriorating progressively until becoming unsuitable as the size of the critical nucleus approaches the interfacial layer thickness. In principle, the data can be used to find more general fitting procedures for the surface tension as a function of bubble curvature see, e.g. 36 for a more complete discussion based on experimental data.
The DI model is highly flexible and, in principle, can be extended to include the effect of dissolved gases in the metastable liquid or heterogenous nucleation over different kinds of surfaces 14 . An important feature of the DI model that was up to now exploited only for model fluids (Van der Waals EoS and Lennard-Jones fluids) is the possibility to describe the dynamic evolution of the cavitation bubble coupled with the inertial effect associated with liquid motion, to the point that even extreme events like shock waves launched in the liquid at bubble collapse, both on free space and near boundaries, can be addressed [37][38][39] . The dynamic approach can be extended in the spirit of Landau and Lifshitz's fluctuating hydrodynamics to include the effect of thermal fluctuations 13,14 and may also account for the wettability of solid boundaries and the coupling with macroscopic flow motion. Finally, the resulting system of stochastic partial differential equations allows to deal with systematic heat injection for heat transfer problem.

Methods
Diffuse interface model. The liquid-vapor system is described via a Diffuse Interface (DI) modelling based on the square gradient approximation (SGA) of the (Helmholtz) free energy functional, proposed in a seminal work by van der Waals 40 . This approach is understood as a particular case of the more general Density Functional Theory (DFT) 33 and allows to control the thermodynamic behavior of the fluid through the choice of effective equations of state expressed in terms of the bulk free energy density, f b (ρ, T) . The total (Helmholtz) free energy functional is where (T) is the temperature dependent capillary coefficient that enables to tune the surface tension of the flat liquid-vapor interface σ 0 (T) of the specific fluid at any given temperature, the surface tension of water 20 in the present work. More precisely, σ 0 (T) is evaluated as: www.nature.com/scientificreports/ expressed both in terms of the equilibrium density profile ρ eq (x) across the flat interface, or directly through the EoS (see 39 for details). In the above equations, ω b (ρ, T) = f b (ρ, T) − µ b (ρ, T)ρ is the grand potential density, µ b (ρ, T) = ∂f b /∂ρ is the chemical potential and ρ sat V /L (T) are the temperature dependent vapor and liquid densities at saturation, respectively. The specific procedure to adapt the IAPWS-95 EoS (the actual expression of free energy f b (ρ, T) used to reproduce the water thermodynamic properties) to the DI approach is described in the next subsection. Free energy barrier evaluation. As discussed in the main text, nucleation is a rare event which calls for specialized algorithm to determine the free energy barrier and, more generally, the minimum free energy path. The total free energy functional is minimized, with the constraint of given total mass, to obtain the equilibrium condition where the external chemical potential, µ ext , and the system temperature are assigned. The above condition implies that the generalized chemical potential, including the capillary term, is spatially homogeneous at equilibrium. When the external chemical potential corresponds to a metastable liquid condition, one possible (unstable) equilibrium solution is the critical bubble immersed in the metastable liquid, described in terms of the critical density profile, ρ c (r) with r the radial coordinate originating at the bubble center. In this work, the powerful string method 24,25 is deployed to obtain the density profile of the critical bubble, corresponding to the saddle point of the (Landau) free energy landscape, i.e. the transition state 13,16 . In a nutshell, the string method numerically approximates the minimum energy path (MEP) describing the continuous sequence of density configurations, ρ(r, α) , along the transition path parametrized by the transition coordinate α . The images ρ k (r) forming the string are evolved over the pseudo-time t according to the steepest-descent algorithm representing a relaxation evolution of Eq. (4), written in spherical coordinates. A reparametrization procedure redistributes the images after evolving Eq. (5) over a single pseudo-time step t . This two-steps procedure is iterated up to the complete convergence of the whole string to the MEP. The configuration ρ c (r) laying on the saddle point corresponds to the critical bubble. The relaxation dynamics of the string images, Eq. (5), is implemented by a centered, second order accurate finite difference scheme, with (pseudo-)time advancement performed by a forward Euler scheme. The forward scheme is used for its simplicity and stability properties, considering the time accuracy is not an issue here since we only need the steady, fully relaxed solution.
The critical radius, the energy barrier and the curvature dependent surface tension, are then measured following 41 as respectively (the temperature T, being a parameter, has been avoided to ease the notation). 19 is here used to accurately describe the thermodynamic properties of water, both in the liquid and the vapor phases. The EoS provides an empirical expression for f b (ρ, T) obtained by fitting a large experimental dataset of stable liquid and vapor states. The values at metastable conditions are extrapolated up to the spinodal states, where the condition ∂p/∂ρ = 0 is met. In the whole unstable region, ρ sp V < ρ < ρ spL , the EoS exhibits spurious oscillations, preventing the direct application into the DI model. In fact, the presence of density intervals where ∂p/∂ρ > 0 in the unstable region, produces unphysical stable-states. The original f b (ρ, T) is corrected by modifying the expression in the unstable region at each temperature:

Equation of state. The IAPWS-95 EoS
is obtained by enforcing the following six conditions requiring the continuity of the pressure and its derivative and of the free energy at the spinodal states: The following third order polynomial supplemented with an exponential function is used as a prototype function for the modified pressure, p mod = Aρ 3 + Bρ 2 + Cρ + D + E exp(ρ)ρ 2 . The cubic expression guarantees that no spurious oscillations occur in the unstable region. The sixth free parameter, F, appears as the integration constant when integrating p mod to obtain the free energy, as f mod b = (p mod /ρ 2 )dρ + F . Other possible choices for f mod b have been tested, showing a very weak sensitiveness on the obtained cavitation pressure ( Supplementary  Information (SI)).

Tolman fit.
An iterative best fitting procedure has been exploited to estimate the Tolman length from the DI surface tension σ * in critical conditions, Eq. (8). For each temperature T, the data consist of the set {σ * i = σ * (µ i lev , T)} , measured at the i-th metastability level ( i = 1, M ). They are plotted, for several temperatures, as a function of the corresponding critical radius, R * i , in Fig. 4. Iterating over n, at fixed T, the Tolman length δ n is obtained by using the Tolman law as fitting function σ * i (T) = σ 0 (T)/(1 + 2δ n (T)/R * i ) , with σ 0 (T) the temperature dependent water surface tension of the flat interface provided by 20 . At the n-th iteration level, the Tolman length δ n is obtained by employing the subset {σ * i (T)} with i from n + 1 to M thus progressively eliminating from the fit the data at the smallest radii. The "corrected"-CNT, with the curvature dependent surface tension, has then been applied to estimate the nucleation rates J n (i) = J n (µ i lev , T) at the different metastable liquid conditions, by employing the different δ n fitting values. The optimal Tolman length, δ opt (T) , is the one that minimizes the L2 norm of the difference between the "corrected"-CNT estimate and the rate measured from the DI model, || log 10 J n − log 10 J DI || = ( M i=1 | log 10 J n (i) − log 10 J DI (i)| 2 ) 1/2 . The fitting error for the different n and several temperatures are reported in the SI. The caption of Fig. 4 reports the optimal Tolman lengths at the different temperatures considered in the plot.
Nucleation rate. The evaluation of the nucleation rate in the different thermodynamic conditions follows the line of the Kramers' theory 32 , providing the transition rate of the nucleation process as the escape of a random walker from the free-energy potential well, k = [( ∪ exp(−β��(x))dx)( ∩ exp(β��(x))/D(x)dx)] −1 , along the transition coordinate x. In this expression β = 1/(k B T) , D(x) is the diffusion coefficient, and ∪, ∩ are shortcuts to indicate that the integration is performed over the the well and barrier basin, respectively. The corresponding nucleation rate is then evaluated as J = Ŵ 0 exp(−�� * /(k B T)) , with the prefactor Ŵ 0 = c 0 k B Tσ 3 /(η�p) derived in 11 in the case of bubble nucleation. The normalization constant c 0 is not uniquely defined in the framework of CNT, but we adopted the expression c 0 = ρ L ρ V /m 227 which compared favorably with the value obtained with MD in 11 . It is worth stressing that the specific expression, formally resembling that of plain-CNT, contains the distinctive feature that the energy barrier, the surface tension and the pressure jump, are all evaluated from the critical density profile ρ c (r) obtained with the string method (Eqs. [6][7][8]. When comparing the rates in Fig. 5, the plain-CNT prediction is obtained by substituting in the above expression the surface tension of the flat interface σ 0 , �� * = 16πσ 3 0 /(3�p 2 ) , and p = p sat − p L . Conversely, in the Tolman-corrected-CNT, the curvature dependence of the surface tension σ (R) = σ 0 /(1 + 2δ/R) is considered in the derivation of the specific expressions 11 , while the value of the Tolman length δ is extracted from our data, as explained in the previous subsection.
p mod (ρ spV ) = p(ρ spV ) p mod (ρ spL ) = p(ρ spL ) ∂p mod ∂ρ (ρ spV ) = 0 License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.