Molecular dynamics simulations in the microcanonical ensemble are performed to study the collapse of a bubble in liquid water using the single-site mW and the four-site TIP4P/2005 water models. To study system size effects, simulations for pure water systems are performed using periodically replicated simulation boxes with linear dimensions, *L*, ranging from 32 to 512 nm with the largest systems containing 8.7 × 10^{6} and 4.5 × 10^{9} molecules for the TIP4P/2005 and mW water models, respectively. The computationally more efficient mW water model allows us to reach converging behavior when the bubble dynamics results are plotted in reduced units, and the limiting behavior can be obtained through linear extrapolation in *L*^{−1}. Qualitative differences are observed between simulations with the mW and TIP4P/2005 water models, but they can be explained by the models’ differences in predicted viscosity and surface tension. Although bubble collapse occurs on time scales of only hundreds of picoseconds, the system sizes used here are sufficiently large to obtain bubble dynamics consistent with the Rayleigh–Plesset equation when using the models’ thermophysical properties as input. For the conditions explored here, extreme heating of the interfacial water molecules near the time of collapse is observed for the larger mW water systems (but the model underpredicts the viscosity), whereas heating is less pronounced for the TIP4P/2005 water systems because its larger viscosity contribution slows the collapse dynamics. The presence of nitrogen within the bubble only starts to affect bubble dynamics near the very end of the initial collapse, leading to an incomplete collapse and strong rebound for the mW water model. Although nitrogen is non-condensable at 300 K, it becomes highly compressed and reaches a liquid-like density near the collapse point. We find that the dissolution of nitrogen is much slower than the movement of the collapsing water front, and the re-expansion of the dense nitrogen droplet gives rise to bubble rebound. The incompatibility of the collapse and dissolution time scales should be considered for continuum-scale modeling of bubble dynamics. We also confirm that the diffusion coefficient for dissolved nitrogen is insensitive to pressure as the liquid transitions from a compressed to a stretched state.

## I. INTRODUCTION

The prediction of multiphase or cavitating flows is inherently complex due to the large span in spatial and temporal scales of the distinct physical phenomena involved, ranging from bubble nucleation at molecular scale to sheet or cloud cavitation at the macroscopic scale.^{1–4} Among them, the collapse of bubbles is of particular interest as it is known to cause erosion in hydroturbines through high-pressure pulses, liquid jets, or local high temperatures generated during collapse.^{5–7} Experiments have led the investigation of bubble collapse. However, the spatiotemporal resolution is restricted by limitations of the imaging instruments,^{8} and the measurements of local extreme conditions (e.g., temperature) within the bubble are often indirect.^{9} In addition to experimental measurements, theoretical analyses and computational studies have contributed to our understanding of bubble collapse. The Rayleigh–Plesset (RP) equation^{10,11} is an ordinary differential equation describing spherical cavity dynamics in an incompressible infinite medium. It is a simple yet powerful tool for deriving many basic results regarding bubble dynamics.^{1,12} Meanwhile, computational fluid dynamics (CFD) simulations can provide more details and are capable of probing more complex systems. However, various assumptions are introduced to account for the multiphase nature of the system, such as the source term describing the vapor production rate in a local volume, leading to different number and/or forms of the transport equations being solved.^{13–15}

For the study of bubble dynamics, molecular dynamics (MD) simulations directly apply classical mechanics to generate trajectories of molecular (particle) systems where the particle interactions are computed using computationally efficient molecular-mechanics force fields. The time and length scales approachable by MD simulations are typically orders of magnitude smaller than those probed by continuum-scale simulations, making MD only suitable for probing molecular scale phenomena, such as nucleation.^{16–19} However, with increasing computing power, recent MD simulations are capable of probing much larger system sizes.^{20–23} In this way, system size effects can be carefully estimated and MD simulation results can serve as a validation for the underlying physical models utilized in continuum-scale simulations. Previous MD simulations have looked into bubble collapse in both Lennard-Jones (LJ) fluids and molecular fluids with or without a nearby surface or external field.^{21–29} To expand the applicability and validate principles seen in MD simulations of aqueous bubble collapse, we study systems with a length scale exceeding a half-micron, compare different water models while considering their thermophysical properties, and examine the effects of non-condensable gas at realistic concentrations.

In this work, we report MD simulation results for the collapse of a bubble in water using both the computationally efficient single-site mW model^{30} and the very accurate four-site TIP4P/2005 water model.^{31} System size effects were studied for both water models, but an order of magnitude larger box length (or three orders of magnitude more molecules) can be accessed with the mW model. Results obtained from MD simulations are compared with numerical solutions from the RP equation, which further explains the different collapse dynamics predicted by the two water models. To study the effects of nitrogen on bubble collapse, a single-site mW-compatible nitrogen model was developed, whereas the TraPPE nitrogen model^{32} without adjustments was used in conjunction with the TIP4P/2005 water model. In addition to system properties, spatially localized properties, such as density and kinetic temperature, were analyzed and visualized with heat maps.

## II. METHODS

### A. Molecular models

#### 1. Single-site models

^{30}utilizes both two-body and three-body (distance–distance–angle) interaction terms such that

*r*

_{ij},

*r*

_{ik}, and

*θ*

_{ijk}are the distance between particles

*i*and

*j*, the distance between particles

*i*and

*k*, and the angle between particles

*j*,

*i*, and

*k*. The mW model uses the following parameters:

*A*= 7.049 556 277,

*B*= 0.6 022 245 584,

*m*= 4,

*n*= 0,

*γ*= 1.2,

*a*= 1.8,

*θ*

_{0}= 109.47°,

*λ*= 23.15,

*ϵ*/

*k*

_{B}= 3114.5 K (where

*k*

_{B}is the Boltzmann constant), and

*σ*= 2.3925 Å. Both two-body and three-body interactions are truncated at a distance of

*aσ*(4.3065 Å).

_{2}–N

_{2}and water–N

_{2}interactions are described by a shifted-force LJ 12-6 potential,

*ɛ*

_{ij}and

*σ*

_{ij}are the well depth and diameter for the unshifted LJ potential. For both N

_{2}–N

_{2}and water–N

_{2}interactions, the truncation distance,

*R*

_{c}, is set to 9 Å. Following prior work,

^{33}fluid phase equilibrium calculations are used to determine suitable LJ parameters. Here, $\u03f5N2\u2212N2$ and $\sigma N2\u2212N2$ are fitted to the critical temperature and density of neat N

_{2}. The resulting parameters are $\u03f5N2\u2212N2/kB=135$ K and $\sigma N2\u2212N2=3.62$ Å, and the model also yields satisfactory saturated vapor pressures (see Fig. S1 in the supplementary material). Given the simplicity of the single-site water and nitrogen models, suitable parameters for the unlike interactions cannot be determined from combining rules. Thus, the $\u03f5water\u2212N2$ and $\sigma water\u2212N2$ parameters are fitted to the solubility of nitrogen in liquid water at two state points where we pick the lowest pressure at

*T*= 298 K and the highest pressure at 343 K from the data set in the work of Chapoy

*et al.*

^{34}Using data at elevated pressure reduces the statistical uncertainty for the phase equilibrium simulations (due to the higher mole fraction of nitrogen in the aqueous phase) and should be more relevant as the bubble collapse leads to compression of the nitrogen-rich region. The resulting parameters are $\u03f5water\u2212N2/kB=123$ K and $\sigma water\u2212N2=3.22$ Å. Details of the fitting procedure for the N

_{2}–N

_{2}and water–N

_{2}interaction parameters are provided in the supplementary material.

#### 2. Multisite models

^{31}and the three-site TraPPE nitrogen models

^{32}are nonpolarizable and use a rigid molecular geometry. The pairwise additive interactions are described by a combination of LJ and Coulomb potentials,

*q*

_{i}and

*ɛ*

_{0}are the partial charge on site

*i*and the relative permittivity of vacuum, respectively. The Lorentz–Berthelot combining rules are used to determine the LJ parameters for the unlike interactions. The Coulomb interactions are computed using the particle–particle particle-mesh Ewald summation technique

^{35}with a relative force error of 10

^{−5}. Both the LJ interaction and the real-space contribution from the Ewald summation are truncated at a distance of 14 Å. To test the performance of the combination of the TIP4P/2005 and TraPPE nitrogen models, the solubility of

*N*

_{2}in water is computed at the same two state points as for the single-site models. The TIP4P/2005 + TraPPE combination underestimates the experimental solubility by a factor of about 1.3 (see Table S1), which is sufficiently accurate for the present study to not warrant any adjustment of the unlike interaction parameters.

### B. MD simulations

The MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package.^{36} However, minor modifications were required for the largest systems containing more than 5 × 10^{8} molecules to change some variables to 8-byte integer, increase the array dimensions, and initialize velocities from the Maxwell–Boltzmann distribution with zero linear momentum for subregions of the system. All simulations were carried out using cubic boxes with periodic boundary conditions. For the setup of an initially two-phase system, the space is divided into a surrounding liquid region and a spherical bubble region with radius *R*_{w}. For the ease of analysis and visualization, the origin of the bubble region is placed at the center of the cubic box. Due to the very low saturated vapor pressures predicted by the mW and TIP4P/2005 water models at 298 K (which underpredict the experimental vapor pressure by factors of 60 and 4, respectively),^{37} the estimated number of water molecules within the largest bubble would be less than 5. Therefore, for the neat water simulations, all molecules are placed inside the liquid region and the bubble is essentially a void. The setup procedure before starting the bubble collapse simulation can be divided into two steps: (1) First, the water molecules in the liquid region are equilibrated using a simulation in the *NVT* ensemble at *T* = 298 K where a spherical LJ wall at *R*_{w} prevents the molecules from entering the bubble region. (2) Next, the LJ wall is removed and the system undergoes a very short relaxation for 0.1 ps in the *NVT* ensemble at *T* = 298 K. This short period is sufficient to remove a temperature spike upon removal of the LJ wall. The setup for the water/nitrogen mixture adopts a similar procedure, except that the initial configuration is obtained by combining a previously equilibrated pure water system containing molecules only outside of the bubble region with N_{2} molecules placed randomly inside the LJ wall—that is, a situation where the nitrogen concentrations in the bubble and in the surrounding liquid region are not in equilibrium.

To probe the bubble collapse dynamics unperturbed from thermostating interferences, following the 0.1 ps pre-equilibration, the simulations are carried out in the microcanonical (*NVE*) ensemble and the timer for the bubble collapse starts at the beginning of the *NVE* simulations (*t* = 0). For simulations with the single-site models, the time step used for equilibration and the majority of the collapse period is 5 fs, but it is reduced to 2.5 fs as the collapse point is approached. Otherwise, the large velocities of particles near the bubble wall (see Sec. III C) would lead to problems with energy conservation. The switch to the shorter time step was made at about 0.9 of the bubble collapse time, *t*_{c}, for the systems started with the larger initial compression and at about 0.97 *t*_{c} for those with the smaller initial compression. For the simulations with multisite models, the time step is set to 1 fs throughout the entire trajectory. Given the very large system sizes used here, the trajectories need to be broken into smaller parts to fit within the wall-clock limits of the computer systems used. To ensure that the trajectories are not perturbed by these restarts, we used binary restart files, kept the grid structure for the Message Passing Interface (MPI) processes (e.g., a 32 × 32 × 32 MPI grid for the mW-512-0.075 system), and retained LAAMPS settings (e.g., fix, neigh_modify, compute, etc.) in all the input files. It was also checked that no warnings were issued in the log file.

For all the systems studied in this work, the overall water density was fixed at values of 998.0 and 995.0 kg/m^{3} for the mW and TIP4P/2005 models, respectively. These values are slightly above (by 0.03% and 0.06%, respectively) the predicted values for the saturated liquid density of these models at *T* = 298 K;^{37} that is, given sufficient time, these systems would relax into a homogeneous liquid phase. Confining the water molecules to the region outside of the bubble leads to an initial state with a higher density in this region and, consequently, higher (system) pressure. For the mW water model, the effect of initial liquid density was explored by performing two series of simulations with *R*_{w}/*L* = 0.15 and 0.075 (where *L* is the edge length of the cubic simulation box), yielding bubble volume fractions of 0.0141 and 0.001 77, respectively, at setup. The corresponding initial pressure values for the systems at *t* = 0 ps are about 800 and 110 atm, respectively. For the TIP4P/2005 water model, simulations were performed only for *R*_{w}/*L* = 0.15, and the initial pressure is about 370 atm. For pure water systems, simulations were performed using various *L* values to study system size effects. For convenience, pure water systems in this work are denoted by a combination of water model, *L* in units of nm, and *R*_{w}/*L*. Simulations for the water/nitrogen mixture were performed for mW-128-0.15 and TIP4P/2005-64-0.15 systems and are denoted as mW-128-0.15-N_{2} and TIP4P/2005-64-0.15-N_{2}, respectively. The nitrogen mole fraction $xN2$ is about 1.2 × 10^{−5}. This concentration is close to the solubility of N_{2} at *T* = 298 K and *p* = 1 atm, and the resulting initial nitrogen pressure inside the bubble region estimated using the ideal gas equation of state is ∼1 atm. As is true for neat water systems, the equilibrium state for the mixture systems is a homogenous liquid phase with all of the nitrogen molecules dissolved—that is, in the condensed state with small partial molar volume. A summary of the systems studied here is provided in Table I.

. | L
. | . | . | . | R_{0}
. | $dRdt0$ . | t_{c}
. | $dRdtBCP$ . |
---|---|---|---|---|---|---|---|---|

Compounds (models) . | nm . | R_{w}/L
. | N_{water}
. | $NN2$ . | nm . | m/s . | ps . | m/s . |

Water (mW) | 32 | 0.15 | 1 093 177 | 0 | 4.75 | −31 | 18.75 | −990 |

64 | 0.15 | 8 745 414 | 0 | 9.56 | −33 | 33.125 | −3120 | |

64 | 0.075 | 8 745 414 | 0 | 4.81 | −7 | 34.875 | −790 | |

128 | 0.15 | 69 963 308 | 0 | 19.16 | −25 | 63.25 | −4060 | |

128 | 0.075 | 69 963 308 | 0 | 9.64 | −11 | 69.625 | −2610 | |

256 | 0.15 | 559 706 465 | 0 | 38.37 | −28 | 123.625 | −8010 | |

256 | 0.075 | 559 706 465 | 0 | 19.26 | −7 | 148.975 | −2840 | |

512 | 0.075 | 4 477 651 718 | 0 | 38.48 | −6 | 309.00 | −4660 | |

Water (mW)/nitrogen (single-site) | 128 | 0.15 | 69 963 308 | 832 | 19.16 | −25 | 63.625 | −1520 |

Water (TIP4P/2005) | 32 | 0.15 | 1 089 891 | 0 | 4.84 | −25 | 69.4 | ⋯ |

48 | 0.15 | 3 678 381 | 0 | 7.21 | −35 | 98.5 | −320 | |

48 | 0.15 | 3 678 381 | 0 | 7.21 | −34 | 98.5 | −400 | |

64 | 0.15 | 8 719 125 | 0 | 9.74 | −3 | 121.5 | −410 | |

Water (TIP4P/2005)/nitrogen (TraPPE) | 64 | 0.15 | 8 719 125 | 104 | 9.71 | −9 | ⋯ | ⋯ |

. | L
. | . | . | . | R_{0}
. | $dRdt0$ . | t_{c}
. | $dRdtBCP$ . |
---|---|---|---|---|---|---|---|---|

Compounds (models) . | nm . | R_{w}/L
. | N_{water}
. | $NN2$ . | nm . | m/s . | ps . | m/s . |

Water (mW) | 32 | 0.15 | 1 093 177 | 0 | 4.75 | −31 | 18.75 | −990 |

64 | 0.15 | 8 745 414 | 0 | 9.56 | −33 | 33.125 | −3120 | |

64 | 0.075 | 8 745 414 | 0 | 4.81 | −7 | 34.875 | −790 | |

128 | 0.15 | 69 963 308 | 0 | 19.16 | −25 | 63.25 | −4060 | |

128 | 0.075 | 69 963 308 | 0 | 9.64 | −11 | 69.625 | −2610 | |

256 | 0.15 | 559 706 465 | 0 | 38.37 | −28 | 123.625 | −8010 | |

256 | 0.075 | 559 706 465 | 0 | 19.26 | −7 | 148.975 | −2840 | |

512 | 0.075 | 4 477 651 718 | 0 | 38.48 | −6 | 309.00 | −4660 | |

Water (mW)/nitrogen (single-site) | 128 | 0.15 | 69 963 308 | 832 | 19.16 | −25 | 63.625 | −1520 |

Water (TIP4P/2005) | 32 | 0.15 | 1 089 891 | 0 | 4.84 | −25 | 69.4 | ⋯ |

48 | 0.15 | 3 678 381 | 0 | 7.21 | −35 | 98.5 | −320 | |

48 | 0.15 | 3 678 381 | 0 | 7.21 | −34 | 98.5 | −400 | |

64 | 0.15 | 8 719 125 | 0 | 9.74 | −3 | 121.5 | −410 | |

Water (TIP4P/2005)/nitrogen (TraPPE) | 64 | 0.15 | 8 719 125 | 104 | 9.71 | −9 | ⋯ | ⋯ |

During the large-scale MD simulations, we compute only the virial pressure of the system from the diagonal components of the system pressure tensor because this does not add to the computational expense. That is, the pressure calculation does not account for spatial variations of the virial in these inhomogeneous systems. Previous works investigating droplets^{38} and bubbles^{39} with relatively small systems have calculated the normal and transverse components of contours as a function of radial distance using the Irving–Kirkwood^{40} and Harasima^{41} pressure tensors. Here, we have not followed this more complex approach because of (a) the additional computational workload, (b) complications arising from the three-body term used for the mW model, and (c) problems of statistical nature due to rapid evolution of the interface during collapse and the fraction of interfacial molecules being negligible for our large systems. The pressure inside of the bubble is estimated here using the instantaneous number densities of (vapor-like) water and nitrogen molecules within the bubble radius, the average kinetic temperature of these particles, and the ideal gas law.

In addition to the MD simulations probing bubble collapse, which are not run until a given system reaches its equilibrium state due to the prohibitive computational expense for these very large system sizes, we also carried out simulations for smaller systems to investigate nitrogen dissolution and evaporation in bubbly systems and nitrogen diffusion in homogeneous systems.

To assess the time scales associated with nitrogen dissolution from an overcompressed bubble and nitrogen evaporation into a bubble from a supersaturated solution, simulations for the single-site models were conducted under conditions where a bubble is thermodynamically stable for the chosen system size, and the N_{2} molecules were initialized at locations distributed either wholly inside or outside of the preformed bubble. The system was comprised of 133 155 mW water molecules and 133 mW nitrogen molecules (a N_{2} mole fraction of 0.001) placed into a cubic box of *L* = 16 nm (yielding a total mass density of 974 kg/m^{3}). The simulations were carried out in the canonical ensemble with the temperature set to 298 K. Based on prior simulations for neat water,^{37} a bubble radius of *R*_{w} ≈ 2.3 nm (*R*_{w}/*L* = 0.14) is estimated to be stable at this condition, and we assumed that the presence of nitrogen leads to a higher pressure inside the bubble and, hence, a larger bubble. The systems were initialized with a preformed bubble of *R*_{w} = 2.4 nm devoid of any water molecules and alternately containing either all 133 N_{2} molecules for the simulations probing dissolution or no N_{2} molecules for the simulations probing evaporation. The remaining water and N_{2} molecules were homogeneously dispersed throughout the remainder of the box volume. A short energy minimization was performed to a force tolerance of 10^{−6} kcal/mol/Å, after which a 25 ps pre-equilibration was conducted with a time step of 0.5 fs and a spherical LJ wall at *R*_{w} in place to prevent water and nitrogen molecules from entering or exiting the bubble region. At this point, the wall was removed and the system was allowed to come to chemical equilibrium over the span of 100 ns. A larger time step of 5 fs was used during this period. The temperature was maintained using the Nosé–Hoover thermostat with a damping time of 100 time steps, and the linear center of mass drift was removed every 100 time steps. Five independent replicas were used for both dissolution and evaporation studies.

_{2}molecule. The cubic box edge length was adjusted to obtain specific densities, using the total mass of water and nitrogen, ranging from 995 to 1005 kg/m

^{3}for the single-site models and from 992 to 1000 kg/m

^{3}for the multisite models. For the chosen system size with

*L*≈ 6.7 nm, the lowest densities are much higher than the maximum density allowed for a stable bubble, i.e., the systems are in the homogeneously stretched phase.

^{37}The temperature of 298 K was maintained using the Nosé–Hoover thermostat with a damping time of 100 time steps. To verify that the 10 N

_{2}molecules will not aggregate together at the higher concentration, test simulations were initialized with all N

_{2}molecules in a single cluster. Within less than a nanosecond, this aggregate dispersed entirely, emphasizing that N

_{2}motions are uncorrelated. As will be seen later, statistically indistinguishable diffusion coefficients were obtained for the 10:10 000 and 1:10 000 N

_{2}:H

_{2}O ratios. Subsequent simulations for the systems with 10 N

_{2}molecules were initialized placing them in dispersed locations. For the single-site models, the initial structure was relaxed for 25 ps with a time step of 0.5 fs, which was followed by a 250 ps equilibration using a 5 fs time step. The diffusion coefficient was computed over a 250 ns production run with a 5 fs time step. For the multisite models, the initial structure was relaxed for 1 ps with a time step of 0.1 fs, followed by a 100 ps equilibration using a 1 fs time step. The production run spanned 20 ns with a 1 fs time step. A total of 12 independent simulation replicas were conducted for the purpose of uncertainty estimation. The average simulation pressure was computed via the molecular virial

^{42}and reported as the average over replicas, with the uncertainty estimated from the standard error of the mean (SEM). Diffusion coefficients were computed via the Einstein relation,

^{42}

^{,}

*D*, is obtained from the mean squared displacement (MSD) (the quantity in angular brackets). The value of

*D*was fitted to MSD data as a function of time, ensuring that the MSD was in the diffusive regime. Because the uncertainty in MSD is heteroscedastic (increasing with time), a weighted linear fit to MSD data from each simulation replica was performed with weights of 1/

*σ*

^{2}, where

*σ*

^{2}is the variance in the MSD at each

*t*.

^{43}The variance is the square of standard deviation at a time lag

*t*=

*τ*, which is best represented as the SEM of the observations that contributed to the average value at

*τ*(because the weights 1/

*σ*

^{2}pertain to a normal distribution, rather than the 3D

*χ*-squared statistical distribution associated with an MSD in 3D space). The estimates for

*D*from each simulation replica were then used to compute an average, and the uncertainty in

*D*is reported as the SEM obtained from the 12 replicas.

### C. Data analysis

The calculation of bubble volume follows a similar protocol as our previous work on equilibrium properties of bubbly water systems.^{37} The analysis uses a mesh-based algorithm with a resolution of 0.2 × 0.2 × 0.2 nm^{3} to divide space into cells. If a cell contains a liquid-like water molecule (i.e., a water molecule that is surrounded by at least two other water molecules using an O–O distance threshold of 0.33 nm) or its center falls within 0.33 nm of at least two other liquid-like water molecules, then the cell is determined to belong to the liquid region. All non-liquid-like cells are considered as vapor-like. The bubble is defined as the largest contiguous cluster of vapor-like cells, where two cells are considered to be part of the same cluster when they share at least one corner. The bubble volume is calculated from the number of cells belonging to the largest cluster. For the water/nitrogen systems, only the water molecules are considered for sorting cells into liquid-like and vapor-like. To analyze the bubble volume efficiently for systems with large number of molecules, a parallelized neighbor search algorithm with domain decomposition is applied.^{44} For consistency, the bubble collapse point (BCP) is defined as the time when the bubble volume decays below a cutoff of 0.08 nm^{3} for the neat water systems (i.e., ten contiguous meshes and roughly equivalent to the volume of three water molecules) or the largest void in the system moves to a new location and as the first volume minimum for the water/nitrogen mixtures. The bubble radius, *R*, is calculated from the bubble volume assuming a spherical shape, which was found to hold reasonably well even for voids with a radius of about 1 nm.^{37} The interfacial bubble wall velocity $dRdt(t)$ is obtained numerically from *R*(*t*) with a center difference scheme for middle points and one-sided difference scheme for boundary points. These analysis methods are used for both the pure water and the water/nitrogen systems. For the mixtures, the fraction of nitrogen molecules contained within the bubble, $\alpha N2$, is also of interest. As will be shown later, the water/nitrogen interface becomes somewhat diffuse after the bubble collapse, and a unique definition is not possible. Here, a nitrogen molecule is considered to be inside the bubble when the distance between its center of mass and the bubble center is less than *R*(*t*).

The time resolutions near the collapse point for the mW and TIP4P/2005 systems are 0.125 and 0.5 ps, respectively. Given the large system sizes used in the current study, it is not possible to store coordinates during the entire trajectory and much of the analysis is performed on the fly. In cases where the coordinates are needed to obtain spatially resolved images near the collapse point, a short part of the trajectory is re-run using the restart file from a few picoseconds before the collapse point, this time writing the coordinates of all molecules to file.

The initial values (i.e., at the start of the *NVE* part of the trajectory, *t* = 0) for the bubble radius, *R*_{0}, and wall velocity, $dRdt0$ are listed in Table I. Due to the short thermalization period after the LJ wall is removed, *R*_{0} is slightly different from *R*_{w}, and $dRdt0$ assumes a small negative value; that is, the bubble radius is shrinking. The simulations for the mW-*L*-0.075 systems show consistently less negative values for $dRdt0$ than the mW-*L*-0.15 simulations started with a more compressed liquid region.

### D. Predictions using the Rayleigh–Plesset equation

*ν*

_{liq}and

*ρ*

_{liq}represent the kinematic viscosity and density of the surrounding liquid;

*γ*represents the surface tension at the bubble interface; and Δ

*P*(

*t*) is defined as

*P*

_{∞}(

*t*) −

*P*

_{bub}(

*t*), which is the difference between the external pressure infinitely far from the bubble and the pressure within the bubble. Without the viscosity and surface tension terms, the RP equation is named the Rayleigh equation. To map MD simulation conditions into input parameters for solving the RP equation,

*P*

_{∞}is taken to be equal to the system pressure,

*P*

_{sys}, while

*P*

_{bub}is assumed to be 0. Although the latter is a satisfactory approximation for the neat water systems, it clearly would not hold for the water/nitrogen mixtures where we find that the nitrogen density inside the bubble becomes quite large as the BCP is approached. For simulations in the microcanonical ensemble, both

*ρ*

_{liq}and

*P*

_{sys}are time-dependent, where

*ρ*

_{liq}is estimated assuming an infinitely sharp interface. The time evolutions of

*ρ*

_{liq},

*P*

_{sys}, and their ratio (last term of the RP equation) for the mW-256-0.15, mW-512-0.075, and TIP4P/2005-64-0.15 systems are shown in Fig. S3. The trajectories of

*ρ*

_{liq}and

*P*

_{sys}are elongated by repeating the last data point when solving the RP equation. The initial conditions [

*R*

_{0}and $dRdt0$] are taken from Table I. The numerical solutions of the RP and Rayleigh equations are obtained using the MATLAB ode45 package.

^{45,46}

It was previously noted that some of the thermophysical properties predicted with the mW and TIP4P/2005 water models differ significantly from each other and sometimes deviate substantively from experiment.^{37} Since the bubble collapse simulations encounter a range of compressed liquid densities, the viscosities for the homogeneous liquid phase represented by the mW and TIP4P/2005 water models were calculated at *T* = 298 K for the densities spanning 995–1016 kg/m^{3} (see Fig. S4). We find that the change in *ν*_{liq} with respect to *ρ*_{liq} is not significant and, hence, *ν*_{liq} is assumed to be constant. Values of 3 × 10^{−7} and 8 × 10^{−7} m^{2}/s are chosen here for the mW and TIP4P/2005 models, respectively.^{37} Although the water molecules near the bubble interface heat up significantly during the bubble collapse process (see Sec. III D), the system temperature does not change appreciably and assuming a constant value for *ν*_{liq} appears reasonable. To test the sensitivity of the RP solutions to the viscosity term, the value of *ν*_{liq} is varied. Based on prior work,^{37}^{,} *γ* decreases by only around 15% when the interface is changed from planar to higher curvature for a roughly spherical bubble with a radius of around 1 nm. Thus, a good first-order approximation is to assume a constant value of *γ*. To test the sensitivity of the RP predictions to *γ*, the value predicted for the flat interface at *T* = 298 K (*γ* = 0.065 N/m for both the mW and the TIP4P/2005 models)^{19,37} and a lower value of 0.050 N/m (accounting for curvature and local heating^{19,37,47}) are considered here.

## III. RESULTS AND DISCUSSION

### A. Vapor bubble collapse for the mW water model: Effects of system size and outside pressure

The time evolution of bubble volume, *V*, bubble radius *R*, and wall velocity, $dRdt$ for pure water systems with the interactions described by the mW water model are presented in Fig. 1. With increasing system size (*L*, compare data with same line style but different color) and initial compression of the outside liquid region (compare solid and dashed lines), the most negative wall velocity, $dRdtBCP$, captured preceding the BCP, becomes more negative, indicating a more violent collapse (for numerical data, see Table I). The most violent collapse is observed for the mW-256-0.15 system, and the corresponding $dRdtBCP$ is about −8 km/s; this corresponds to Mach 3.5 (based on the value of 2.3 km/s for the sound speed of liquid water at *T* = 298 K and *P* = 1 atm calculated using the mW model and the method of Lustig^{48}). For this system, a clear rebound is observed with the bubble growing back to *R* > 3 nm (*V* > 100 nm^{3}). For the mW-128-0.15 and mW-512-0.075 systems, $|dRdtBCP|$ exceeds 4 km/s, and a minor rebound to *R* ≈ 0.9 nm (*V* ≈ 3 nm^{3}) is observed, but the rebound is very short-lived with the bubble vanishing after 12 ps in both cases.

To assess system size effects, the time evolution of bubble volume and radius is plotted in reduced units in Fig. 2, where bubble volume and radius are normalized by their initial values at *t* = 0, and time is normalized by the corresponding bubble collapse time, *t*_{c}. Note that this type of nondimensionalization process is a common practice in fluid mechanics,^{49} and bubble collapse dynamics solved on larger scales are often presented in this way as well.^{50,51} For both the mW-*L*-0.15 and mW-*L*-0.075 systems, a converging behavior for the initial bubble collapse process is observed as the system size increases. As mentioned above, the rebound after initial collapse is only observed for a few systems but, where rebounding is observed (i.e., the mW-128-0.15 and mW-256-0.15 systems), there is remarkable consistency between system sizes in reduced units. It should be noted that, under these conditions, the RP equation is not capable of predicting a rebound after the initial collapse point.

To estimate the limiting behavior (*L* = *∞*) for the bubble collapse dynamics, an extrapolation based on inverse linear dimension was performed. First, the reduced time range (0 ≤ *t*/*t*_{c} ≤ 1) was divided into intervals of length 0.01. For each time slice, a linear fit was performed for *V*/*V*_{0} as a function of 1/*L*, and the *y*-intercepts yield the data for *L* = *∞* (see Fig. 3). The correlation coefficients are satisfactory with *R*^{2} > 0.995 for 0.1 ≤ *t*/*t*_{c} ≤ 0.8. Somewhat larger deviations from linearity are observed at the beginning and near the BCP, which may be explained by the slightly different initial conditions in terms of $dRdt0$ (see Table I) and the extreme heating and supersonic wall speeds observed for only the larger systems. The estimates of reduced bubble volume and radius for *L* = *∞* are also shown in Fig. 2. For both the mW-*L*-0.15 and mW-*L*-0.075 systems, the collapse trajectories for the largest bubble are very close to the limiting behavior, indicating that the collapse behavior of a bubble with initial radius around 40 nm is already a good representative for the collapse of mesoscopic bubbles with radii at least on the order of micrometers, which can be observed experimentally with high-speed cameras.

### B. Differences between water models and reconciliation through the RP equation

The time evolution of the bubble volume, *V*, bubble radius, *R*, and wall velocity, $dRdt$, for pure water systems predicted by the TIP4P/2005 water model are presented in Fig. 4. With the initial liquid pressure value for the TIP4P/2005-*L*-0.15 systems falling in between those for the mW-*L*-0.15 and mW-*L*-0.075 systems (see Fig. S3) and using the same simulation setup, one would expect a similar collapse behavior for the mW and TIP4P/2005 models. However, comparison between Figs. 1 and 4 indicates qualitatively different dynamics. The simulations for the TIP4P/2005 model predict a local minimum in $dRdt$ at the early stage of bubble collapse process. Near the BCP, a local maximum, though less obvious, is also observed in $dRdt$. These two local extrema correspond to two inflection points in *R*(*t*) where the sign of curvature changes, whereas the mW water model predicts a continuously decreasing $dRdt$ (increasing in magnitude) before reaching the BCP. Because of the elongated collapse process, the bubble collapse observed for the TIP4P/2005 water model is overall much less violent.

To understand the stark difference between these water models, we turn to the RP equation and recognize model-dependent thermophysical properties.^{37} By evaluating the signs of different terms on the left-hand side of the RP equation, we observe that it is possible to reach an inflection point in *R*(*t*), where $d2Rdt2=0$, when the negative viscosity term (due to $dRdt<0$ for a collapse process) balances out the other terms that are positive (due to *ν*_{liq}, *γ*, and Δ*P* being positive). As noted earlier, the mW and TIP4P/2005 models yield kinematic viscosity values at ambient conditions of 3 × 10^{−7} and 8 × 10^{−7} m^{2}/s, respectively (see also Fig. S4). The corresponding experimental value is 8.9 × 10^{−7} m^{2}/s.^{52} Thus, the TIP4P/2005 model underestimates the viscosity by about 10%, whereas the mW model underestimates it by a factor of 3.

Figure 5 presents MD simulation results and numerical solutions of the Rayleigh and RP equation for the largest systems: mW-256-0.15, mW-512-0.075, and TIP4P/2005-64-0.15. Although the Rayleigh equation captures the initial bubble dynamics (i.e., when Δ*P* is very large) reasonably well for these three systems, it fails to predict the qualitative differences observed for the mW and TIP4P/2005 water models that emerge later in the collapse process. Inclusion of the surface tension term leads to a faster collapse, whereas an opposite effect is observed for the viscosity term. This physically makes sense as reducing the surface area is an additional driving force, while viscosity describes the resistance to flow. Overall, the agreement between MD simulation results (shown in black) and the RP equation utilizing the appropriate kinematic viscosity and surface tension for a given water model (shown in cyan) is very satisfactory, which has also been reported in some previous studies.^{25,28,53}

For the mW-256-0.15 system that yields the most violent collapse, the Rayleigh equation with the additional surface tension term yields slightly better agreement with the MD simulations than the Rayleigh or RP equations, but the differences are rather small when the appropriate *ν*_{liq} for the mW model is used. When *ν*_{liq} = 10^{−6} m^{2}/s (i.e., close to the experimental value) is used, then an inflection point is also observed in the RP prediction for the mW-256-0.15 system, but this qualitatively fails to match the MD trajectory. For the mW-512-0.075 system, the RP equation with the appropriate *ν*_{liq} for the mW model yields the best agreement with the MD trajectory, and the Rayleigh equation with the surface tension term slightly underpredicts the time until the BCP. Again, using *ν*_{liq} close to the experimental value yields an inflection point and extended collapse process that is not in agreement with the MD results. Reducing the value of surface tension to 0.05 N/m while keeping *ν*_{liq} = 3 × 10^{7} m^{2}/s leads to slightly smaller wall velocities and a slight increase in *t*_{c}.

For the TIP4P/2005-64-0.15 system, only the RP equation with appropriate *ν*_{liq} for the water model predicts collapse dynamics with an inflection point in *R*(*t*) that matches the MD trajectory. That is, a larger *ν*_{liq} value is required to yield the inflection point and the stretched collapse dynamics.^{25} Using a smaller value for *γ* results in additional stretching of the collapse trajectory. In contrast, the Rayleigh equation and the RP equation with *ν*_{liq} = 2 × 10^{−7} m^{2}/s yield a much more violent collapse process. However, the RP equation with the appropriate *ν*_{liq} value is not capable of predicting a second inflection point near the BCP that leads to an acceleration of the final segment of collapse. Prior simulations in the canonical ensemble indicate that changes in surface tension^{19,37} and viscosity^{37} are relatively small even for (metastable) vapor bubbles with a radius of about 2 nm. We surmise that the acceleration of the final part (*R* < 2 nm) of the collapse process may be caused by thermodynamic forces as bubbles are not even metastable anymore.^{18,37}

Bubble dynamics obtained for the TIP4P/2005-*L*-0.15 systems were also examined in reduced units (see Fig. S5) to assess whether the system sizes studied here are sufficient to obtain the limiting behavior (*L* = *∞*). Linear fits were made for the data at the three system sizes studied for the TIP4P/2005 water model (see Fig. S6). Unlike the findings for the mW systems, where a very satisfactory linear correlation was observed for most of the bubble collapse process, the linear fits for the TIP4P/2005-*L*-0.15 systems are good (*R*^{2} > 0.95) only at intermediate times (0.33 ≤ *t*/*t*_{c} ≤ 0.83), and *R*^{2} < 0.5 is found for *t*/*t*_{c} < 0.27 and for *t*/*t*_{c} in the vicinity of 0.92. To test the sensitivity to initial conditions, a second independent run was performed for the TIP4P/2005-48-0.15 system, but both runs yielded very similar data (see Fig. S5). It appears that the two *t*/*t*_{c} regions with the worst *R*^{2} values are near the two inflection points occurring at somewhat different reduced times for these TIP4P/2005-*L*-0.15 systems. It is possible that the linear correlation breaks down for bubble collapse processes with inflection points in *R*(*t*) or that larger system sizes are needed to reach the regime where a linear correlation holds. With the currently available computing resources, it is not feasible to run simulations for *L* = 128 and 256 nm for the TIP4P/2005 model, which has more interaction sites per molecule and requires a larger cutoff and Ewald summation together with a shorter time step, while at the same time exhibiting much slower collapse dynamics due to its larger (but closer to experiment) viscosity.

### C. Temporally and spatially localized extreme conditions

To quantify the temperature and density distributions in time and space, the system is divided into spherical shells containing ∼1000 water molecules. Here, only the translational kinetic temperature is used to allow for a fair comparison of the simulations employing the single-site and multisite models. At the beginning of the collapse process, the temperature is near 300 K for all shells, and the density profile shows strong oscillations reflecting the presence of the wall (see Fig. S7). As the bubble shrinks and the wall velocity increases, the radial temperature profile develops an exponentially decaying behavior with the innermost shell being always found near the highest temperature for any of the shells (see Fig. S7). The change in density profiles is less pronounced (apart, of course, from changes in the location of the bubble wall), but there is a trend toward a slightly increased density ($\u22481020$ kg/m^{3}) near the bubble wall at *t*/*t*_{c} = 0.95. This transient density enhancement may also affect the surface tension.

We focus here on the innermost shell because these molecules are most likely to cause material damage, cavitation erosion, and sonoluminescence. The time evolutions of the average temperature of the innermost shell (*T*_{1000}) and the threshold temperature for the ten molecules with the highest translational kinetic temperature (*T*_{10}) for the pure water systems are presented in Fig. 6. In general, the time evolutions of *T*_{1000} and *T*_{10} are qualitatively similar for the same system as one would expect for a Maxwell–Boltzmann distribution, and *T*_{10} ≈ 4 *T*_{1000} holds reasonably well. For the mW water systems, a sharp peak is observed in both *T*_{1000} and *T*_{10} near the collapse point. With increasing system size and initial outside liquid pressure (comparing mW-*L*′-0.15 and mW-2 *L*′-0.075 systems starting with the same initial bubble radius), a higher peak temperature in *T*_{1000} and *T*_{10} is observed. A comparison of Figs. 1 and 6 leads to the conclusion that a more violent collapse [i.e., larger value of $dRdtBCP$] corresponds to a stronger heating effect for the pure water systems. However, the local heating is more pronounced. As the system size is doubled, we observe that $dRdtBCP$ also approximately doubles (see Table I), whereas *T*_{1000} quadruples. For the largest bubble size studied in this work, *T*_{1000} peaks at values exceeding 20 000 and 6000 K for the mW-256-0.15 and mW-512-0.075 systems, respectively, which is about 70 and 20 times higher than the initial value of *T*_{1000} ≈ 300 K (see Fig. S7). These extreme temperatures well in excess of 1000 K agree with experimental observations.^{5,54} Our simulation data support the Jarman hypothesis^{55} that sonoluminescence is thermal in origin, which is based on the empirical premise that sonoluminescent flashes are broad spectrum and would require emitters with temperatures of several thousand Kelvin.^{56}

Compared to the mW systems, a qualitatively different heating behavior is observed for the TIP4P/2005 systems studied here. First, the system size does not seem to have an effect on the evolution in reduced time for both *T*_{1000} and *T*_{10}. Second, both *T*_{1000} and *T*_{10} increase slightly at the beginning of the collapse process, but then their time evolution flattens out at *t*/*t*_{c} ≈ 0.2. At the transition point, the corresponding absolute times are around 14, 20, and 24 ps for *R*_{0} = 4.8, 7.2, and 9.6 nm, respectively. These times correspond approximately to the initial inflection point in *R*(*t*) (see Fig. 4). In other words, after passing the first inflection point when viscosity effects begin to dominate and $dRdt$ starts to increase, the local heating of the innermost shell is no longer observed.

Figure 7 presents the space–time evolution of the specific density and kinetic temperature in the vicinity of the BCP for the mW-128-0.15 and TIP4P/2005-64-0.15 water systems. Within the 2 ps preceding the collapse point, the bubble radius decreases by around 4 nm for the mW-128-0.15 system, while a decrease of less than 1 nm is observed for the TIP4P/2005-64-0.15 system, again reflecting a much more violent collapse for the mW water system with the lower viscosity. For both water models, the moving interfacial region has an intermediate density and a width of around 0.5 nm. As the collapse point is approached, a slight decrease in the interfacial width is observed. For the mW-128-0.15 system, an extremely high density exceeding 1.6 g/cm^{3} is found at the BCP for shells within 0.6 nm of the center, and a wave of higher density shells moving away from the bubble center is observed for the first 0.5 ps after the BCP. The mW-128-0.15 system exhibits a small rebound (see Fig. 1), but the interface of the rebounding bubble is rather diffuse. The space–time evolution for the kinetic energy shows that the heating is most pronounced at the interface with a well-defined hot spot with temperatures as high as 9000 K at the BCP for the mW-128-0.15 system. Interestingly, the temperature does not show a feature similar to the radiating density wave immediately following the collapse. Cross-sectional heat maps near the BCP are shown for this system in Fig. S8. Before the BCP, the highest temperatures and voxel velocities (all pointing inward) are observed for the interfacial region. After the BCP, however, the region of highest temperature remains at the center, whereas the aftermath of the collapse leads to high density and voxel velocity (all pointing outward) waves radiating away from the center. The collapse process for the TIP4P/2005-64-0.15 system is much less violent, and the space–time evolution shows a more diffuse picture without a large spike in density nor temperature at the BCP (see Fig. 7).

### D. Bubble dynamics in water/nitrogen systems

To investigate the effects of a non-condensable gas on bubble collapse dynamics, we also performed simulations for water/nitrogen mixtures using single-site and multisite models. It is important to emphasize that these simulations were initialized with all of the inert gas molecules placed inside the bubble but that the number of nitrogen molecules is close to the solubility limit and the initial pressure in the bubble is only about 1.2 bars. Figure 8 presents the bubble dynamics results for the water/nitrogen systems and a comparison to the corresponding pure water systems. In the presence of nitrogen, complete bubble collapse is not observed for the mW-128-0.15-N_{2} and TIP4P/2005-64-0.15-N_{2} systems, and the minimum bubble radii are found to be 1.5 and 1.2 nm, respectively. Interestingly, the effect of nitrogen on the bubble dynamics is negligible for the overwhelming majority of the trajectory. For the mW-128-0.15-N_{2} system, the time of 63.6 ps until the first minimum in bubble radius is only 1% longer than the *t*_{c} for the mW-128-0.15 system. Even for the TIP4P/2005-64-0.15 systems, the bubble radii and wall velocities closely trace each other for about 110 ps (i.e., about 0.9*t*_{c} for the system without nitrogen). The close agreement between neat water systems and water/nitrogen mixtures for the majority of the collapse trajectory also extends to the density in the liquid region, the system pressure, and their ratio (see Fig. S9) that constitutes the last term on the right-hand side of the RP equation. We also note that the bubbles maintain spherical symmetry throughout the initial collapse process (see Fig. S9).

However, the effect of nitrogen on the final part of the initial collapse process and the subsequent dynamics is very pronounced. For the mW-128-0.15-N_{2} system, the terminal wall velocity is nearly three times smaller in magnitude than for the corresponding pure water system (see Table I). For the former system, the initial collapse is followed by a much stronger and longer lasting rebound (see Fig. 8). Eventually, the rebounding bubble shrinks again, followed by another rebound, but the second and third collapses are much softer. The *R*(*t*) values for the second and third minima are 2.8 and 2.6 nm, respectively. Due to computer time limitations, the trajectory for this large system (containing 7 × 10^{7} water molecules) was stopped after 200 ps because it became obvious that reaching an equilibrium state with no bubble and all nitrogen molecules dissolved and well distributed in the liquid would be unattainable. For the TIP4P/2005-64-0.15-N_{2}, the presence of nitrogen leads to a dramatic dampening of the initial collapse process; that is, there is not a clear minimum in *R*(*t*) and also not a second inflection point in $dRdt$. We observe a trend toward a decreasing bubble size for the last 50 ps of the trajectory for this system but, again, approach toward the equilibrium state with no bubble appears to be extremely slow compared to the initial collapse process (see inset of Fig. 8). Prior large-scale simulations presented in the work of Shekhar *et al.*^{21} also investigated the effects of a non-condensable gas, but the initial gas (number) density in the bubble was more than three orders of magnitude larger than in the present work, i.e., the equilibrium state presented in the work of Shekhar *et al.* study would be a two-phase system.

To provide insight into the pronounced effects on the final phase of the bubble collapse caused by the presence of nitrogen, we calculated a few properties for only the nitrogen molecules contained within the bubble (see Fig. 9). For the mW-128-0.15-N_{2} system, the fraction of nitrogen molecules within the bubble, $\alpha N2$, decays only slowly from near unity (831 out of 832 nitrogen molecules contained within the bubble at *t*_{0}) to a value of 0.9 during the first 50 ps of the trajectory despite that the bubble volume shrinks by a factor of 6. This is followed by a steep decline to a value of 0.43 when the first minimum is reached; that is, there is some degree of mixing during the more violent part of the collapse process. However, the bubble rebound is accompanied by demixing, and $\alpha N2$ reaches 0.84 before the bubble starts to shrink again. Subsequent bubble size oscillations are also accompanied by mixing as the bubble volume decreases and demixing as the bubble volume increases. Thus, the mixing appears to be mostly limited to an interfacial zone (see below).

The density and temperature for nitrogen molecules in the bubble also undergo oscillations driven by initial collapse and subsequent rebounds. For the mW-128-0.15-N_{2} system, the specific density of nitrogen inside the bubble, $\rho N2$, reaches in excess of 1000 kg/m^{3}, thereby matching the specific density of the surrounding liquid; however, accounting for the difference in molecular mass, the number density of nitrogen is smaller than that of water. During the first rebound, $\rho N2$ falls to about 20 kg/m^{3}, but quite surprisingly it only falls to 200 kg/m^{3} during the second rebound. The reason is that the first rebound reaches a bubble volume that is about ten times larger than that for the second rebound, whereas $\alpha N2$ decays by only 10%. The temperature of nitrogen molecules within the bubble, $TN2$, reaches maxima of 6400, 570, and 450 K during the first, second, and third bubble collapses, respectively, and falls to 300 K during the rebounds. Previous hybrid RP-MD studies also predicted extreme heating exceeding 10 000 K near the collapse point for violent bubble collapse with inert gas content.^{57,58} Using the density and temperature values, we can also estimate the pressure in the bubble. For convenience, Fig. 9 presents data using the ideal gas equation of state. [We have performed constant-*NVT* MD simulations for pure nitrogen systems at selected $(\rho N2,TN2)$ values observed during the initial collapse for the mW-128-0.15-N_{2} system and find that the ideal gas equation of state holds well for all but the most extreme state points, but it underestimates the pressure by up to a factor of 4 for $TN2>3000$ K. Since nitrogen’s Boyle temperature is $\u2248325$ K,^{59} the ideal gas equation of state holds rather well for $TN2<500$ K.] The simulations indicate that, at the BCP of the mW-128-0.15-N_{2} system, nitrogen is compressed to a density comparable to that observed for a pressure in excess of 1 GPa and pressure equivalents of 370 and 480 bars for the second and third collapses, respectively.

As discussed before, the collapse process for the TIP4P/2005-64-0.15-N_{2} system is much less violent without a well-defined BCP. This also leads to qualitatively different trends for the properties of nitrogen within the bubble (see Fig. 9). The fraction of nitrogen molecules contained within the bubble decreases at a near constant rate for the first 130 ps of the trajectory and, when the bubble volume becomes steady, $\alpha N2$ fluctuates around 0.7. Since the majority of nitrogen molecules are still contained in the bubble, but the volume has decreased by a factor of about 600, a very high nitrogen density ($\rho N2\u2248500$ kg/m^{3}) and pressure ($pN2\u2248450$ bars) are also reached for the TIP4P/2005-64-0.15-N_{2} system. Despite the less violent collapse dynamics, it is still much faster than the time required for the dissolution of nitrogen molecules. As also seen for the corresponding pure water system, the temperature of nitrogen molecules only increases until the first inflection point in *R*(*t*), and $TN2$ does not exceed 400 K.

Figure 10 presents the space–time evolution of water density, nitrogen density, and temperature near the initial collapse point for the mixture systems. The V-shaped interface observed in the water density heat map for the mW-128-0.15-N_{2} system indicates a clear rebound, consistent with the bubble dynamics results presented in Fig. 8. Preceding the collapse, the interfacial width remains at around 0.5 nm, while it grows rapidly during the rebound. Cross-sectional heat maps near the BCP for the mW-128-0.15-N_{2} system are shown in Fig. S10. As expected from the soft collapse for the TIP4P/2005-64-0.15-N_{2} system, the radius of the bubble and the interfacial width barely change within the ±2 ps near the first minimum in *R* at *t* = 134 ps. For both the single-site and the multisite models, an incomplete mixing is observed near the collapse point, where water molecules do not penetrate deeply into the interior of the bubble and where nitrogen molecules are highly compressed within the bubble. Interestingly, for the mW-128-0.15-N_{2} system, the radius of the region of high nitrogen density remains nearly constant after the collapse, while the water front is receding. Instead the nitrogen density in this dense core is decreasing. These observations reflect a difference in the time scales between the slow dissolution of nitrogen into the surrounding liquid water and the fast collapsing water front (see Sec. III E). With respect to temperature, the presence of nitrogen leads to slightly more extreme temperatures for the mW-128-0.15-N_{2} system than the mW-128-0.15 system, while no significant difference is observed for the TIP4P/2005 systems.

### E. Nitrogen dissolution, evaporation, and diffusion

As discussed in Sec. III D, the simulations indicate that the bubble collapse process for the water/nitrogen mixture systems does not lead to destruction of the bubble and dissolution of the nitrogen molecules; instead we observe a long-lived bubble containing a very highly compressed nitrogen gas. The simulations for the bubble collapse used a very large system size (see Table I) that does not allow us to follow the trajectories for more than a few hundred ps. Thus, we conducted simulations for the single-site models using a 8^{3} smaller system (in regard to number of molecules and volume) to assess the time scales that are pertinent to reaching an equilibrium state from initial states representing an overcompressed bubble or an oversaturated solution region. For the system studied here, the initial state with all N_{2} molecules placed inside the bubble corresponds to a pressure of about 100 bars inside the bubble and the initial state with all N_{2} molecules dissolved in the liquid region corresponds to a N_{2} mole fraction of 10^{−3}.

The trajectories for the systems characterized by mostly dissolution or evaporation are shown in Figs. 11 and S11, respectively. During the first 200 ps of the dissolution process, the bubble rapidly grows to *R*_{w} = 2.44 nm and loses about 20% of the N_{2} molecules. That is, the gas bubble with the highly compressed N_{2} gas prefers a slightly larger volume than a neat water system with the same number of water molecules. This phase of the process reflects the transition from a perfectly spherical bubble with a sharp interface between liquid and gas found for the initial configuration (and enforced by a retaining wall at *R* = 2.4 nm) to a slightly less spherically symmetric bubble with an interfacial width similar to the molecular size (about 0.3 nm) and all of the N_{2} molecules that are nominally outside of the bubble still being close to the interface. The radial density profile calculated over the first 1 ns shows some N_{2} enrichment at the bubble side of the interface (see Fig. 11).

Similarly for the evaporation process simulations, the bubble radius rapidly shrinks over the first 200 ps to *R*_{w} = 2.22 nm while less than two N_{2} molecules nominally enter the bubble (see Fig. S11). That is, the oversaturated liquid phase requires a larger volume (i.e., smaller *R*_{w}) than the neat water system with the same number of water molecules. The radial density profile calculated over the first 1 ns does not exhibit N_{2} enrichment at the bubble the interface (see Fig. S11).

After the initial rapid changes, the bubble is seen to shrink or grow more slowly over time as N_{2} molecules dissolve into or evaporate from the liquid region, respectively. After about 40 ns, the dissolution and evaporation trajectories reach the same equilibrium state with *R*_{w} = 2.33 nm and $\alpha N2=0.36$. This equilibrium state corresponds to a pressure inside the bubble (calculated from the N_{2} density using the ideal gas law) of about 40 bars and a solution mole fraction of about 6 × 10^{−4} that are in agreement with the bulk partitioning experiments^{34} (see also Fig. S2).

The radial density profiles indicate that the slight enhancement of N_{2} number density at the bubble side of the interface is also present when equilibrium is reached. For the dissolution process, the N_{2} density in the liquid region shows a decrease with increasing *r* at 5 and 10 ns, but the liquid region profiles flatten by 20 ns. For the evaporation process, the N_{2} density profiles in the liquid region are reasonably flat throughout the trajectory. Based on the density profiles in the liquid region not showing a strong gradient well before the bubble reaches equilibrium, we surmise that diffusion is not rate-limiting for the dissolution/evaporation processes.

Figure S12 shows the evolution of system pressure, the pressure inside the bubble, and the surface tension. During the dissolution process of the compressed gas inside the bubble, the system pressure and the pressure in the bubble decrease, and the reverse holds for the evaporation process from the supersaturated solution. Remarkably, the surface tension does not show pronounced changes during the dissolution and evaporation processes, indicating that the surface tension is not very sensitive to the presence of nitrogen over the range of nitrogen densities within the bubble observed during these trajectories, and the limiting value reached at *t* = 80 ns falls within 1% of the value obtained from equilibrium simulations for vapor bubbles in neat water.^{37}

As should be expected, the individual trajectories do not follow exactly the same path, but their dissolution and evaporation processes are quite similar. The rapid oscillations that are seen for both *R*_{w} and $\alpha N2$ are largely due to thermal fluctuations, the resolution of the mesh-based bubble detection algorithm, the conversion of mesh volume to *R*_{w} assuming a perfect spherical shape, and partitioning of the system into liquid and bubble regions without allowing for an interfacial region. A simple exponential fit of the form *A* exp(−*t*/*τ*) + *C* to *R*_{w}(*t*) describes the evolution of the average values of *R*_{w} and $\alpha N2$ remarkably well and reveals a characteristic time constant *τ* of 1.0 × 10^{1} ns for both evaporation and dissolution (see Figs. 11 and S11).

That is, for the systems investigated here, the characteristic time scale for dissolution is two orders of magnitude slower than the time scale for the initial bubble collapse. This is consistent with the observation that there is insufficient time for N_{2} to equilibrate with the surrounding liquid during the bubble collapse, giving rise to a rebounding effect for the more violent collapse observed for the mW-128-0.15-N_{2} system, a softening for the collapse of the TIP4P/2005-64-0.15-N_{2} system, and a long-lived bubble dissolving only on much longer time scales for both systems. The nonequilibrium rates of bubble growth/shrinkage due to dissolution and evaporation from/into a bubble are found here to be comparable, suggesting that traversing the liquid/vapor bubble interface and mass transport away from/toward the interface in the solution phase are rate-limiting. Of course, the extent of separation of the collapse and dissolution time scales may depend on the boundary conditions governing the collapse, but we surmise that it will play a pivotal role during the more violent phase of gas bubble collapse.

As mentioned above, gas diffusion in the solution phase contributes to the separation of time scales between bubble collapse and gas dissolution. Experimental measurements of the diffusion coefficient of dissolved gases in water are challenging because of the low solubility of common gases. In fact, a popular approach to determine gas diffusion coefficients is by observing the collapse dynamics of large bubbles (with diameters larger than 100 *μ*m) usually at ambient pressure.^{60} The Taylor dispersion method uses a laminar flow tube, and detection of the dispersion curve of the sparingly soluble gases is aided by performing the measurements at elevated pressure.^{61} Neither of these methods allows for measurements at pressures approaching the saturated vapor pressure of water. This has led to some uncertainty regarding the pressure dependence of the diffusion of dissolved gases, and some researchers studying bubble growth have used the Chapman–Enskog relation^{62} to extrapolate diffusion coefficients into the low-pressure regime,^{63} while others have assumed constant diffusion coefficients independent of pressure.^{64} For MD simulations, there are no conceptual problems with studying the tracer diffusion in a homogeneous solution phase, but the small number of gas molecules in the simulation box leads to a statistical challenge.

Data for N_{2} diffusion in water at 298 K covering a range of system densities are presented in Fig. 12. The N_{2} diffusion coefficients obtained for the multisite models exhibit no discernible density dependence while the pressure ranges from −16 to 160 bars, i.e., from the homogeneously stretched to compressed liquid regime. This finding is in agreement with the experimental measurements using the Taylor dispersion method presented in the work of Cadogan *et al.*^{61} that yield *D* = (1.99 ± 0.05) × 10^{−5} cm^{2}/s for pressures from 120 to 450 bars. Averaging the N_{2} data from the simulations for the multisite models yields *D* = (1.97 ± 0.08) × 10^{−5} cm^{2}/s. Here, it should be noted that the TIP4P/2005 water model also yields an excellent diffusion coefficient at ambient conditions.^{31,37} Prior experimental measurements of the N_{2} diffusion coefficient give values of *D* = (3.1 ± 0.3) × 10^{−5} cm^{2}/s (interpolated from data at 293 and 303 K) using the bubble collapse method,^{60}^{,} *D* = (2.30 ± 0.12) × 10^{−5} cm^{2}/s (interpolated from data at 293 and 303 K) using a stagnant liquid layer method,^{65} and *D* = (2.01 ± 0.20) × 10^{−5} cm^{2}/s using the dispersion technique.^{66}

The simulation data for the single-site models yield somewhat lower statistical uncertainties and indicate that the N_{2} diffusion coefficient may decrease by slightly less than 2% as the pressure increases from −155 to 365 bars, but compared to the data presented in the work of Cadogan *et al.*^{61} and multisite model data, the N_{2} diffusion coefficient of *D* = (3.59 ± 0.02) × 10^{−5} cm^{2}/s at |*p*| ≤ 5 bars is overestimated substantially by a factor of 1.8 (as one may expect from the over prediction of the self-diffusion coefficient in the neat water phase by a factor of 2.9 for the mW water model^{30,37}). Based on the experimental data presented in the work of Cadogan *et al.*^{61} and the present simulation study, it is clear that the Chapman–Enskog relation should not be applied to describe the pressure dependence of the diffusion of a dissolved gas. This is because the mean free path length for gas in the solution phase is dictated by the solvent density, which barely changes as a function of pressure for a liquid solvent far from its critical point. The closely scattered data points near zero pressure found here for the single- and multisite models demonstrate that the diffusion coefficient of the dissolved gas does not diverge as the liquid transitions from a compressed to a stretched state. Thus, the observation that the experimentally measured rate of bubble growth upon reduction in pressure can be satisfactorily described under the assumption that the diffusion coefficient is inversely correlated with pressure^{63} may instead indicate that the growth is not entirely diffusion-controlled.

As also observed in our prior study for neat water,^{37} the simulations for the mixture indicate a statistically significant change of 1% or 2% for the water diffusion coefficients described by the multi- or single-site models, respectively, over the range of densities investigated here (see Fig. S13). The decrease in water diffusion rate with decreasing pressure can be attributed to the more ice-like behavior exhibited at lower densities.^{37}

## IV. CONCLUSIONS

Molecular dynamics simulations are used to probe the collapse dynamics for pure water and water/nitrogen mixtures. For the computationally efficient mW water models, simulations were performed for a range of system sizes (*L* ranging from 32 to 512 nm) and two values of the initial bubble volume fraction (*α* = 0.0141 and 0.001 77), yielding initial values for the system pressure of 800 and 100 bars, respectively. We find that, under these conditions, an initial bubble radius of 40 nm is already a good representative for the collapse of micrometer-sized bubbles. Due to the underprediction of viscosity by the mW model, the collapse process is rather violent with a strong rebound for the larger systems and temperatures in excess of 5000 K for the 1000 water molecules belonging to the innermost shell (located at the liquid/vapor interface) of the collapsing bubble. The simulations for the TIP4P/2005 water model, which accurately predicts the water viscosity, yield much elongated collapse dynamics where inflection points corresponding to local extrema in $dRdt(t)$ and only modest heating until the initial inflection point are observed. To understand the difference, MD simulation results are compared with the numerical solutions of the RP equation, and the overall agreement is very satisfactory when all terms in the RP equations are based on the properties of the model systems. It is found that inclusion of the viscosity term in the RP equation is necessary for observing the initial inflection point in *R*(*t*), at which the viscosity effect starts to dominate.

The most important finding of this study pertains to the striking differences observed for water/nitrogen mixtures compared to the neat water systems. Even under conditions where all gas molecules would be dissolved in the homogeneous solution phase once equilibrium is reached, we find a clear separation of the time scales for bubble collapse and gas dissolution. Here, the initial collapse is incomplete and results in a gas droplet that is compressed to liquid-like densities. The re-expansion of the dense droplet gives rise to a strong rebound followed by additional incomplete collapse/soft rebound cycles. Based on simulations for smaller systems, we surmise that the time scale for gas dissolution is two orders of magnitude slower than the time scale for the collapse process. We recommend that gas partitioning [i.e., the process of converting a non-condensable gas into a dissolved (condensed) species] and an explicit rate for dissolution/evaporation processes be included in computational fluid dynamics simulations investigating two-phase flow. As an aside, we show that the diffusion coefficient of dissolved gas molecules is insensitive to pressure as the solution phase transitions from a compressed liquid (i.e., large positive pressure) to a homogeneously stretched liquid (i.e., large negative pressure) state.

## SUPPLEMENTARY MATERIAL

The supplementary material for this article includes the following: (1) a PDF file containing a description of the development of the single-site nitrogen model, the nitrogen solubility values in aqueous solution using the multisite models, plots of the input data for the Rayleigh–Plesset equation, plots probing the system size effects for the TIP4P/2005 systems, plot of the radial temperature profile at different t/tc values, cross-sectional heat maps near the collapse point for the mW-128-0.15 and mW-128-0.15-N2 systems, plots for the nitrogen evaporation dynamics from oversaturated solution, and the diffusion of water in water/nitrogen mixtures and (2) two MP4 files that show cross-sectional heatmaps for the water and (if present) nitrogen density for systems (a) mW-128-0.15 and mW-128-0.15-N2 and (b) mW-256-0.15.

## ACKNOWLEDGMENTS

This work was supported by the United States Office of Naval Research under Grant No. ONR N00014-17-1-2676 with Dr. Ki-Han Kim and Dr. Yin-Lu Young as the program managers. C.K. was supported by the Office of Science, U.S. Department of Energy, under Contract No. DE-AC02-06CH11357. Part of the computer resources was provided by the Minnesota Supercomputing Institute and the Argonne Leadership Computing Facility. We thank Ryan DeFever and Edward Maginn for cross-validation of the shifted-force LJ potential and Joe Katz for alerting us to the controversy regarding the dissolved gas diffusion rate at low pressure.

The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan. http://energy.gov/downloads/doe-public-access-plan.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jingyi L. Chen**: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). **Jesse L. Prelesnik**: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). **Buyun Liang**: Investigation (supporting). **Yangzesheng Sun**: Investigation (supporting); Methodology (supporting). **Mrugank Bhatt**: Conceptualization (supporting); Formal analysis (supporting); Writing – review & editing (supporting). **Christopher Knight**: Software (supporting); Writing – review & editing (supporting). **Krishnan Mahesh**: Conceptualization (equal); Funding acquisition (equal); Writing – review & editing (supporting). **J. Ilja Siepmann**: Conceptualization (lead); Funding acquisition (equal); Methodology (equal); Supervision (lead); Visualization (supporting); Writing – review & editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*The Micro-World Observed by Ultra High-Speed Cameras*

*Fluid Dynamics of Cavitation and Cavitating Turbopumps*

*Computer Simulation Using Particles*

*Computer Simulation of Liquids*

*NVT*simulation of the isobaric heat capacity, speed of sound and Joule–Thomson coefficient

*Fluid Mechanics*

*CRC Handbook of Chemistry and Physics*

_{2}and N

_{2}in water at temperatures between 298.15 K and 423.15 K at pressures up to 45 MPa