Numerical Frameworks for Laser-Induced Cavitation: Is Interface Supersaturation a Plausible Primary Nucleation Mechanism?

Crystallization has been observed in laser-induced cavities in saturated solutions, but the mechanisms behind nucleation of crystals are not entirely clear. A hypothesis is that high solution super...


INTRODUCTION
Crystallization is the formation of crystalline solids from a fluid phase, and the process is utilized as an essential separation and purification technique in, for example, the pharmaceutical and fine chemicals industry. 1 During the process of crystallization, it is necessary to control nucleation of the new crystalline phase in order to produce crystals with the desired quality features, such as crystal size, morphology, or purity. Such a control is far from being trivial, and various technologies have been proposed in an effort to achieve the required crystal quality. 2−4 A promising technology for good spatiotemporal control of the nucleation process is the laser-induced nucleation (LIN) method. 5 By applying a laser pulse to a supersaturated solution, nucleation has been observed in experiments using non-photochemical laser-induced nucleation (NPLIN) 6 and those experiments that induce nucleation by the formation of a clearly observable vapor bubble, or cavity. 7,8 The mechanisms behind the nucleation are not entirely clear for any of the LIN methods, but it is hypothesized that solvent evaporation may play a key role in the nucleation process. 9,10 In the NPLIN method, the solvent evaporation may occur around laser-heated nanoparticles, and in the laser-induced cavitation method the solvent is evaporated at the vapor bubble interface. To investigate the hypothesis that solvent evaporation is a prerequisite for the nucleation, we focus in this work on the laser-induced cavitation method where the vapor bubbles have been observed in experiments and the solvent evaporation is known to take place.
By sufficient laser irradiation of a liquid with a high spectral absorbance, the liquid becomes superheated, and a vapor bubble is formed. The bubble grows explosively due to evaporation of the superheated liquid, and the phenomenon is termed thermocavitation. 11 In liquids with low absorption coefficients, cavitation can be achieved if the laser irradiation reaches the threshold for optical breakdown of the liquid. During optical breakdown, nonlinear light absorption due to a cascading ionization process produces a hot, rapidly expanding plasma that can reach pressures greater than 7 GPa. 12 This phenomenon is called optical cavitation. 13 For crystallization purposes, it is favorable to suppress the required laser energy and the resulting solution temperature due to possible thermal degradation of the solute. 14,15 In this work, we focus on crystallization by laser-induced thermocavitation due to the potential of the latter for nucleation at the temperatures lower than those in optical cavitation, where hot plasma is produced. Thermocavitation also has the possibility for a precise control of the local solution temperature. 11,16 The hypothesis about the phenomena leading to the observed nucleation is that the evaporation of a solvent at the liquid−vapor interface causes an increase of the concentration of solute and, simultaneously, cools the solution. 9,10,17 Both phenomena lead to an increased supersaturation of the solution. This increase may be a prerequisite for primary nucleation since nucleation is more probable with increased supersaturation. 4 To investigate if this hypothesis is plausible, the solution supersaturation needs to be estimated during the laser-induced thermocavitation event. Because of very fast dynamics and small scales of the problem, it is not trivial to experimentally measure that property.
Numerical simulations of the laser-induced thermocavitation event can resolve the solution properties and give estimates of the degree of supersaturation in the solution around the bubble. By comparing such simulations with the bubble dynamics and crystallization observations from experiments, the nucleation hypothesis can be tested with a good degree of certainty. To perform such a test we propose a methodology that, as a first step, involves laser-induced thermocavitation experiments at increasing laser pulse energies. With the increasing pulse energy, the resulting thermocavitation bubble will grow more rapidly and become larger due to the faster evaporation of the solvent. Above a threshold energy, the evaporation should, in a suitable solution, produce an interface supersaturation high enough for nucleation to occur. From these experiments, it is possible to determine the threshold laser pulse energy that produces thermocavitation with observable crystallization, similar to the studies by Yoshikawa et al. 18,19 The second step of the methodology involves numerical simulations of the thermocavitation bubbles produced in the experiments. With an agreement of the bubble radius evolution in the simulations and the experiments, the former can accurately predict the interface supersaturation. If the predicted supersaturation reaches significant levels near the threshold laser energy for crystallization in the experiments, it seems very probable that the hypothesis about high interface supersaturation due to evaporation is indeed the phenomenon that triggers the nucleation.
Several experimental studies on laser-induced nucleation have been conducted, but they typically induce the vapor bubble by optical cavitation. 15,18,20−22 The optical cavitation, together with hot plasma formation and extreme dynamics, are highly nonlinear processes and not trivial to simulate. To carry out the proposed methodology, experiments are needed with laser-induced thermocavitation at increasing laser pulse energies and the resulting event of crystallization. To the best of our knowledge, no such experimental study has been performed. Still, we can test the plausibility of the nucleation hypothesis by using simulations and comparing the obtained results with the data obtained in a case that includes thermocavitation and crystallization at a single laser pulse energy.
A suitable experimental study with the laser-induced thermocavitation and the observed crystallization is the one carried out by Soare et al. 10 The authors used an aqueous solution of an inorganic salt, (NH 4 ) 2 SO 4 , where ink was added to enhance the laser irradiation absorption. The solution was placed between two glass slides, 50 μm apart, and the laser irradiation locally heated the solution to an estimated temperature of 494 K. 23 The bubble diameter rapidly exceeded the distance between the glass slides, whereupon the bubble grew in a seemingly two-dimensional manner in the plane of the glass. About 1 s after the laser pulse, crystals became visible in a ring around the laser focal point. In this ring, optical disturbances were observed after only 30 μs that were suggested to arise by the nuclei that had grown large enough to become visible. In the work of Soare et al., 10 the length scales of the nucleation process were too small to be visually observed and, therefore, the nucleation probably took place before the optical disturbances occurred. The authors assumed that the nucleation took place at the bubble interface during the maximum rate of evaporation. This maximum occurred at the early bubble growth phase. Therefore, we argue that the simulations of this early growth phase are essential to capture the conditions that may lead to the nucleation.
Simulations of similar types of bubbles have been reported in the literature. Magaletti et al. 24 and Zein et al. 25 studied collapsing nanobubbles with simulations that included the effects of phase change and predicted bubble radius evolutions in good agreement with a theoretical model and experiments. With a cavitation model, Sagar and el Moctar 26 studied the collapse phase of a laser-induced vapor bubble, with the predicted bubble shape and interface dynamics in fair agreement with experiments. The cavitation model in that study did not include inertia and surface tension, the phenomena that are of high relevance when the growth phase of a small vapor bubble is to be studied. Koch et al. 27 performed simulations of laser-induced bubbles and took into account the effects of inertia and compressibility, but not those of phase change. Using the initial amount of gas in the bubble as a fitting parameter, their predicted bubble radius evolution was in good agreement with experiments for both the growth and collapse phases of the bubbles. Soare 23 developed a simplified model to describe their laser-induced nucleation experiments. The model assumed four well-mixed fluid regions with uniform process conditions. Heat and mass transfer between the regions were considered, and nucleation and crystal growth were modeled in the liquid region closest to the vapor bubble. The modeled bubble radius evolution, number of crystals and crystal size were in good agreement with the experiments. The model was, however, dependent on fitting parameters, and it is not certain that the model and those specific parameters could be used to model other laser-induced nucleation cases.
In this work, we tested the proposed nucleation hypothesis with the following course of action: first, we performed a multiphase DNS simulation of a laser-induced thermocavitation bubble. We resolved all relevant physical phenomena in the process, such as inertia, compressibility, surface tension, and phase change. For that purpose, we used the numerical framework presented in our previous study, 28 but extended it here with an improved formulation of the interfacial energy transfer and with consideration of the solute transport in the liquid around the vapor bubble. This framework is briefly presented in section 2.1, and we used it to perform a 2Daxisymmetric simulation of the laser-induced vapor bubble reported in the experiment by Soare et al. 10,23 The purpose of the simulation is to capture the early growth phase of the vapor bubble and to quantify the degree of supersaturation of the solution around the bubble. The results from this simulation are presented in section 3.1. Our aim is to show that the hypothesis about high interface supersaturation is plausible for the thermocavitation experiment where nucleation is observed.
Second, to obtain a wide parameter space for testing the nucleation hypothesis, we formulate a comprehensive 1D numerical framework that also accounts for relevant fluid dynamics effects, phase change, and the transport of solute. The reason for introducing the latter framework is the enormous computational cost of multiphase DNS simulations (even for 2D cases) due to the great temporal and spatial resolutions required to accurately resolve the small scales and fast dynamics of the problem. We justify the use of the 1D framework with the following reasoning: we know that it is the solvent evaporation that causes the increase of solute concentration around the bubble. By varying the laser pulse energy and the focus of the laser beam, it is possible to control the superheated solution temperature and the volume of the superheated solution. These solution properties govern the rate and duration of the evaporation. In addition to the laser pulse parameters, the diffusivity of the solute also affects the solute concentration, and the solubility of the solute determines the degree of the supersaturation. To gain a deeper understanding of the nucleation mechanisms and how they can be utilized, we need to investigate the effects of the laser pulse parameters and the solution properties on the interface supersaturation. The use of the 1D framework makes it feasible that the values of the relevant parameters are chosen freely and varied independently. Still, we have to note here that if the thermocavitation bubble is affected by surrounding boundaries, simulations in 2D or 3D are needed to capture the correct behavior. Again, the computational cost of such DNS simulations is prohibitive in a way that such simulations are not yet (and will most likely not be in the foreseeable future) a feasible option for a systematic study that includes a larger set of variations within the space of the aforementioned parameters.
In summary, we aim with the 1D framework to show that the nucleation hypothesis related to the presence of high interface supersaturation, due to evaporation, is plausible as a means for nucleation under a range of realistic and industrially relevant conditions. We will examine the effects of the variations of the solute diffusivity, solute solubility, laser pulse energy, and the superheated liquid volume on the solution supersaturation. A great number of simulations are required for such a study, and in practice this can only be done by using a 1D framework. This framework is presented in section 2.2.
We have organized our study that uses the 1D framework into two investigation cases, with the results presented in sections 3.3 and 3.4. In the first one, we look at how the solute diffusivity and solubility affect the solution supersaturation. We aim at presenting guidelines that will make it possible for any future user of our framework to identify if a specific solution attains conditions favorable for primary nucleation, given a set of laser pulse parameters.
In the second investigation case, we examine how the laser pulse energy and the liquid volume that the energy is absorbed in affect the solution supersaturation in a commonly used aqueous ammonium sulfate solution.
We identify the laser pulse parameters that are suitable for high supersaturation of the given solution, and we show the types of laser pulse parameters that are appropriate for supersaturation at lower solution temperatures. The latter is crucial for solutes sensitive to thermal degradation. Finally, we discuss our main findings and draw conclusions in section 4.

METHODOLOGY
In this section, we present our numerical frameworks and describe the simulations that we conduct. The 2D numerical framework is outlined in section 2.1, and in section 2.1.3 we describe our 2D-axisymmetric simulation of the laser-induced thermocavitation bubble. In section 2.2, we present the 1D framework for spherical symmetric thermocavitation bubbles, and in sections 2.2.2 and 2.2.3 we define the two previously mentioned investigation cases.
2.1. 2D Numerical Framework. To test the nucleation hypothesis about high interface supersaturation, we perform a 2Daxisymmetric multiphase DNS simulation of the laser-induced thermocavitation bubble reported in the experiment by Soare et al. 10,23 with observed crystallization. We use a numerical framework that is based on the volume of fluid (VOF) method. 29 The fluid conservation equations of mass, momentum, and energy are solved together with that for a phase indicator function that tracks the relative position of the two phases. To account for mass and energy transfer between the phases, a phase change model is implemented into the VOF framework by a technique developed by Hardt and Wondra 30 and extended by Kunkelmann. 31 In this technique, the mass and energy transfer are described by a distribution of source terms in the conservation equations of mass and energy. These source terms are located in a region close to the bubble interface and ensure that the correct amount of mass and energy is transported between the phases. The rate of the phase change is determined by the model of Schrage. 32 We have validated the numerical framework by comparing the growth of a spherical vapor bubble in an infinite superheated liquid to the analytical solution by Scriven. 33 The validation was presented in our previous work 28 and showed the two growth rates in good agreement. The governing equations and a detailed description of the phase change model can be found in the Supporting Information 1.1 and 1.2, respectively. In our 2D simulation, where the bubble radius reaches almost 100 μm, it is not computationally affordable to resolve the mass boundary layer of about 100 nm. The problem of very thin mass transfer boundary layers in DNS simulations of bubbles has been addressed in different ways. For example, Aboulhasanzadeh 35 implemented a subgrid model for the boundary layer around a rising bubble where the concentration profile was given by a parabolic function. In our case, the concentration at the bubble interface is not constant, but governed by the evaporation rate that varies significantly during the simulation. This makes it difficult to implement an analytical solution for the transient concentration profile. Bothe et al. 36 solved a onedimensional diffusion problem along radial rays, pointing outward and starting at the bubble interface, to compute the mass flux from the bubble surface to the surrounding solution. We adopt a similar approach by considering averaged values of the phase change rate and liquid conditions across the entire bubble interface. The liquid conditions at the interface are taken as the averaged values over the isosurface C = 0.99 where C is the phase indicator function in the VOF method, and C = 1 corresponds to pure liquid. By using averaged values, it is possible to compute a conservative estimate of the solute concentration in the bubble interface normal direction. This reduces the problem to a 1D case that allows for the required spatial resolution. The solute transport equation is solved in a coordinate system that is normal to, and moving with, the bubble interface. The equation is defined as Crystal Growth & Design pubs.acs.org/crystal Article where c is the solute concentration (in this work given as the mass fraction m solute /m solution ), x is the distance from the bubble interface, and u n is the liquid velocity, relative to the bubble interface. The relative velocity is caused by the evaporation and determined by the incompressible continuity equation according to the Supporting Information 1.3. We assume the solute nonvolatile and therefore prescribe null convective and diffusive solute flux boundary conditions at the interface as cu n | x=0 = 0 and | = . In the liquid far away from the interface, the solute concentration is the initial uniform value c(x = ∞) = c ∞ . Equation 1 is solved using the cell-centered finite volume method in a coordinate system that moves with the bubble interface. The computational grid used for solving eq 1 comprises moving and deformable control volumes defined in the Supporting Information 1.3.
2.1.2. Problem Description. The geometry of the experimental thermocavitation event is illustrated in Figure 1. 10 In the experiment, a supersaturated aqueous solution of an inorganic salt, (NH 4 ) 2 SO 4 , was placed between two glass slides. The solution contained magenta ink to increase the absorbance of the laser irradiation. Perpendicular to the glass, a 6 ns laser pulse was applied with a diameter of about 20 μm. The solution absorbed a part of the laser irradiation and became superheated. Then, a vapor bubble formed that grew explosively, reaching a maximum radius of about 700 μm in the plane of the glass. The bubble interface attained a velocity of well above 20 m/s and the total event, from the bubble formation to collapse, took about 200 μs. The solvent evaporation reaches a maximum after the laser pulse when the superheated liquid at the interface is at its highest temperature. All of the superheated liquid is then evaporated during the early bubble growth phase, and this evaporation causes a maximum of solute concentration that is, after the evaporation has stopped, diffused away from the interface. Therefore, we focus in the 2D simulation on the early bubble growth phase where these phenomena may induce high supersaturation levels.
2.1.3. Simulation Setup. The setup of our 2D DNS simulation is somewhat simplified as compared to the experimental case. The geometry and laser pulse duration are kept as in the experiment. The laser pulse energy is chosen to produce the solution temperature of 494 K after the laser pulse that was estimated in the experimental study. 23 Since this temperature is lower than the spinodal limit of water, around 578 K, the vapor bubble probably nucleated heterogeneously on an existing nanobubble or a particle. The laser pulse energy is assumed uniformly absorbed in a cylinder of liquid, with the laser beam diameter, since the absorbance coefficient of the laser irradiation and the laser beam geometry are not known. Measurements of the absorbance coefficient and distribution of the laser pulse energy would be needed to avoid the assumptions about the conditions in the solution after the laser pulse. These assumptions may be the reason for a possible lack of agreement between the simulation and experimental data.
The material property data of the aqueous ammonium sulfate solution are not readily available for the temperatures and concentrations treated in this simulation. Therefore, pure water properties had to be assumed and temperature dependencies are implemented based on the formulations issued by the International Association for the Properties of Water and Steam (IAPWS). 37−40 These properties include the liquid kinematic viscosity μ L , liquid thermal conductivity k L , surface tension coefficient σ, saturation temperature (as a function of pressure) T sat (p), saturation pressure p sat , specific heat capacity c p and the latent heat of phase change .
The surface tension effects are accounted for by the Continuum Surface Force model by Brackbill et al. 41 The vapor phase is assumed an ideal gas, and the liquid phase is considered compressible according to the Tait equation of state. In the phase change model, a value for the accommodation coefficient χ is needed that is not easily determined.
Marek and Straub 42 have reported values in the range of 10 −3 < χ ≤ 1 for water in different experiments. In this work, a value of χ = 1 is chosen, and we will show that it leads to simulation results that are in fair agreement with experimentally observed bubble growth rates. The choice of χ = 1 is also the most challenging value numerically (for stability reasons), which is useful for assessing the robustness of the numerical framework. 30 It is important for the accuracy of the phase change model that the advection scheme of the phase indicator function does not smear the interface. For that reason, the PLIC interface reconstruction scheme is chosen that maintains a sharp interface throughout the simulation. 43 The time step size is variable and determined based on the fluid velocities and the phase change rate according to the Supporting Information 1.5.
Because of axial symmetry, we adopt a 2D axisymmetric computational domain discretized with a Cartesian grid. The domain is illustrated in Figure 2. The height between the r-axis and the glass wall is 25 μm, and the distance between the z-axis and the outlet is 1800 μm. In the part of the domain between the z-axis and the location where r = 400 μm, the computational grid is uniform, with a cell size of 50 nm. In the outer part of the domain, the grid is increasingly coarsened to reduce the computational cost. A total number of 2.23 million cells are used, and a grid convergence study is performed and presented in the Supporting Information 1.4.
A no-slip and a constant-temperature boundary condition is applied at the glass wall and the symmetry boundary conditions at the r-and z-axes. The temperature at the glass wall is 293 K except at the 10 μm closest to the z-axis where the temperature is 393 K due to laser irradiation. 23 At the outlet boundary, a pressure-outlet condition is specified with a constant pressure of 101 300 Pa.
In this simulation, the nucleation of the vapor bubble is not included, and a finite size bubble is initialized with a radius of 1 μm. Initially, the pressure in the bubble is the gauge Laplace pressure of 117 800 Pa, and the temperature is the corresponding saturation temperature of 396 K. The liquid inside the laser beam, r < 10 μm, is uniformly heated by source terms in the energy equation during the first 6 ns of the simulation. During this time, the liquid temperature increases from the initial 293 K to 494 K. The initial solution supersaturation is σ = 0.005, which corresponds to a solute concentration (mass fraction) of c ∞ = 0.433.  A Cartesian structured grid is used with the vapor bubble initialized at the lower left corner. A symmetry boundary condition is applied to the radial and axial axes. At the glass wall, a no-slip wall boundary condition is prescribed, and at the outlet a pressure-outlet boundary condition is adopted that represents conditions far away from the bubble.
Crystal Growth & Design pubs.acs.org/crystal Article Because of the glass walls, the 2D numerical framework is needed to capture the correct bubble dynamics. Inherently, the 2D framework also predicts local variations of the interface solution properties. These features make the framework suitable for detailed studies of specific thermocavitation events. Still, the computational cost of such DNS simulations is rather excessive. For a systematic study of how the solution properties and laser pulse parameters affect the solution supersaturation, a less computationally costly, but still sufficiently accurate, method is needed.
2.2. 1D Numerical Framework. We now shift our attention to developing a 1D numerical framework in order to investigate the proposed nucleation hypothesis under varying solution properties and laser pulse setups. A relatively low computational cost allows us to examine entire relevant parameter spaces of the solution properties and laser pulse setups. We now consider a spherical symmetric laserinduced vapor bubble evolving in an unbounded domain. All material properties are those for pure water, and temperature dependencies are implemented as in section 2.1.3. A 1D framework was used by Akhatov et al. 44 to successfully model the collapse phase of laserinduced vapor bubbles. Here, we use our 1D framework to simulate the entire bubble lifetime. The fluid region is divided into two zones: the compressible gas phase inside the bubble, 0 ≤ r < R, and the incompressible liquid phase outside the bubble, R ≤ r < ∞, where r is the radius from the bubble center and R is the bubble radius.
Accurate predictions of the pressure and temperature on both sides of the interface are essential during the laser pulse and bubble growth period to get accurate evaporation rates. The properties within the vapor phase change rapidly during and after the laser pulse. Therefore, we use a compressible formulation of the conservation equations for the vapor phase. This ensures that we capture the correct dynamics and properties of the vapor during the entire bubble lifetime. The conservation equations of mass momentum and energy for the vapor phase are where the index v denotes the vapor phase, ρ v the vapor density, u v the radial velocity, p v the pressure, components of the shear stress the total energy per unit mass, ϵ v the internal energy and λ v the thermal conductivity. The viscous effects described by τ are often negligible compared to the other terms in eq 3 and 4 but are included here for increasing the numerical stability of the framework when handling such fast dynamics in the simulations. The vapor pressure and internal energy are modeled as a van der Waals gas and determined by where B v = 458.9 J/(kg K), γ = 1.3, b 1 = 1.694 × 10 −3 m 3 /kg, and b 2 = 1708.6 Jm 3 /kg 2 are water vapor parameters. 44 The rate of phase change at the bubble interface is determined by the model of Schrage 45 that uses concepts from the kinetic theory of gases to compute the flux of fluid molecules across the interface from the temperature and pressure of the phases as i k j j j j j j y where j is the mass flux at the bubble interface, χ is the accommodation coefficient discussed in section 2.1.3, T l is the liquid temperature, and Γ v is a correction factor given by The inclusion of the correction factor Γ v is important if the gas phase departs appreciably from its equilibrium condition. 45 Such a departure certainly occurs during the intense heating during the laser pulse. At the bubble interface, the latent heat of phase change is accounted for by the boundary condition where is the specific latent heat, evaluated at the interface temperature. The temperature of the interface is still an unresolved issue, but in this framework we assume a continuous temperature across the interface as 46 The vapor and liquid velocities differ at the interface due to phase change according to The pressure is also discontinuous across the interface due to the surface tension and viscous effects. The jump condition is implemented as where γ is the surface tension coefficient and μ l the dynamic viscosity. At the bubble center, symmetry boundary conditions are specified for all variables.
From the conservation of mass in the incompressible liquid, we have for the spherically symmetric case where we include the effect of the phase change rate j on the liquid velocity. This effect is often neglected for high density liquids, but becomes significant at high evaporation rates attained after the laser pulse. The evolution of the bubble interface is given by the generalized incompressible Rayleigh−Plesset equation 47 (note that here we also include the effects of phase change on the liquid velocity) as where p ∞ is the pressure in the liquid far away from the bubble. This equation is obtained by substituting eq 16 into a momentum equation for the liquid in the radial direction and integrating it from the vapor side of the interface to far away in the liquid. To compute the pressure in the liquid phase p l (r), as a result of the interface motion, a similar Crystal Growth & Design pubs.acs.org/crystal Article approach can be used where the integration is performed from radius r in the liquid phase to far away in the liquid, according to The phase change rate, and thereby the bubble growth rate, is governed by the liquid temperature at the bubble interface. To determine the interface temperature, we need to accurately resolve the temperature profile within the entire liquid region. This is achieved by solving the following form of the energy conservation equation where α l is the liquid thermal diffusivity, P las is the energy source from the laser pulse, and c p,l is the specific heat capacity of the liquid. The boundary condition far away from the bubble is T l | r=∞ = T ∞ , and at the bubble interface the boundary conditions are given by eq 11 and 12.
We model the liquid heating from the laser pulse using the energy source term, P las (r). This source term represents the laser power density that is absorbed by the liquid solution and applied during the laser pulse, 0 ≤ t ≤ t las . The absorption of the laser irradiation is assumed negligible within the vapor phase, so that P las (r) is applied only in the liquid region (R ≤ r < ∞). We further assume that the irradiation profile of the laser beam is Gaussian according to where r las is the radius of the laser beam, and P max is the power density in the center of the beam. The laser power density, P las , is assumed constant during the laser pulse duration, t las , so that the total absorbed laser energy density, e las , can be determined by Here, e max = t las P max is the laser energy density at the center of the beam, and its value is chosen so that a given total laser pulse energy E las is absorbed in the liquid region, around the initial vapor bubble, as where R 0 is the initial bubble radius. Since the bubble radius R(t las ) ≈ R 0 during the short laser pulses we investigate in this work, the initial radius R 0 is used instead of R in the above integration. The solute concentration in the liquid is computed according to with the same boundary conditions as described in section 2.1.1. The vapor phase governing eqs 2−4 is solved in a stationary coordinate system with the origin at the bubble center, while the liquid temperature and solute concentration equations, eqs 19 and 23, are solved in a moving coordinate system with the origin at the bubble interface. All conservation equations are discretized using the cellcentered finite volume method with the implicit second-order Crank− Nicolson scheme in time.
In the thermocavitation events of this study, the bubble radius may change by more than 3 orders of magnitude. To accurately resolve the fluid dynamics during the entire bubble evolution process, an adaptive grid refinement is implemented. The latter ensures that a sufficiently refined grid is used by splitting or joining cells in areas with high or low gradients and curvatures of the fluid properties. With this approach, we reduce the computational cost while preserving accuracy and avoiding unnecessarily strict time-step limitations due to overly refined grids. The time step size is also variable and chosen according to the Supporting Information 2.1.
In all the simulations, the initial temperature and pressure of the liquid phase are T ∞ = 293 K and p ∞ = 101 000 Pa, respectively. The nucleation of the vapor bubble is not considered in this framework, but it most probably occurs either heterogeneously at, for example, an existing nanobubble or particle, or homogeneously if the spinodal limit for water, at around 578 K, is reached. The initial bubble radius is chosen as R 0 = 0.5 μm, and the vapor is initialized with the equilibrium absolute pressure of 392 230 Pa and the corresponding saturation temperature of 416 K.
In the following section, we verify the implementation of the 1D framework and validate the simulation results against experimental data.
2.2.1. Validation of 1D Framework. We here look at the predicted interface solution supersaturation from the 1D and 2D numerical frameworks to ensure that the frameworks predict similar results and are implemented correctly. More information on this part of our study can be found in the Supporting Information 2.3. The 1D framework is validated against two thermocavitation experiments. The first validation case is presented in this section, and the second case is included in the Supporting Information 2.2. In the first validation case, we compare the predicted bubble dynamics with an experimentally observed laser-induced thermocavitation event reported by Quinto-Su et al. 48 In the experiment, the bubble was induced close to a glass wall on one side so that the domain can be considered semi-unbounded. The bubble dynamics is, however, similar to the fully unbounded domain that we assume in the 1D framework. The experimental study used an ink solution, but in the simulation we assume pure water properties, and the duration of the laser pulse is set to that of the experiment, t las = 200 × 10 −15 s. The values of the laser pulse energy and the laser beam radius used when comparing the simulation results with the experimental bubble evolution are E las = 3.9 μJ and r las = 11 μm. The results from the first validation case are given in section 3.2 and show a good agreement with the experimental data.
The agreement of the bubble evolution verifies that realistic evaporation rates are predicted since it is evaporation that governs the bubble growth rate. The accurate evaporation rate also ensures that realistic values of the interface solute concentration and supersaturation can be computed for this thermocavitation event. Therefore, this simulation setup will be used in our first investigation case where we examine the effects of varying the solution properties on the solution supersaturation.
2.2.2. Case 1: Influence of Solution Properties on Supersaturation. In our first investigation case, we use the 1D framework to study how the solute diffusivity and the solubility affect the supersaturation of the solution around a laser-induced thermocavitation bubble. All other material properties of the solution in the simulation are those for pure water so that the solute concentration is treated as a passive scalar and does not influence the bubble dynamics. Therefore, we can investigate the effects of the solute diffusivity and solute solubility independently. The results from this investigation are given in section 3.3.
We use a fixed set of laser pulse parameters, given in section 2.2.1, that give a bubble evolution in agreement with an experimentally observed laser-induced thermocavitation event. The agreement of the bubble evolution makes it possible to compute realistic estimates of the interface solute concentration and supersaturation for a given set of solution properties.
The solubility of a solute can be a highly nonlinear, discontinuous, and multivariable function. In this study, we are not aiming for a comprehensive investigation of all possible solubility behaviors but focus on the most common features. Therefore, we consider the solubility of the solute to be a function of the temperature according to where T is the solution temperature, T 0 = 273.15 K is a constant, and a, b and n are parameters. The initial concentration of the solution, before the laser pulse, is uniform and at saturation conditions. Thus, We want to investigate how the parameters in eq 24 affect the supersaturation in eq 25. To reduce the number of independent parameters, we nondimensionalize the concentration and the liquid temperature and introduce the solubility parameter C 1 as Equation 25 can then be rewritten according to Here, the only parameters are the solubility parameter, C 1 , and the degree of the temperature dependence, n. Three values of n = (0.7, 1.0, 1.3) are investigated to cover a broad range of possible values and to establish qualitative trends in the results. These values represent solutes with a negative, zero, or positive curvature of the solubility curve, respectively. In this study, we consider C 1 ≥ 0, which represents solubility curves that increase with, or are independent of, the temperature. The minimum value, C 1 = 0, inserted into eq 29 results in σ = c* − 1 and, at large enough values, we have C 1 ≫ 1 and , with eq 29 changing to Thus, the limits of σ at the inputs C 1 = 0 and C 1 = ∞ are both independent of C 1 . Because of this independence, it is sufficient to investigate an intermediate part of the total parameter space, C 1 ∈ (0, ∞), that has a significant effect on the relative supersaturation. From this investigation, we show the influence of all shapes of the solubility curve, in the form of eq 24 and for a given n, on the solution supersaturation. To completely define our investigation case of important solution properties, we also need a relevant parameter range for the solute diffusivity.
The evaporation of the solvent increases the solute concentration of the solution and produces a concentration gradient around the bubble. This gradient induces diffusion of the solute away from the interface and the rate of the diffusion mass flux is proportional to the diffusion coefficient D AB . On the basis of the diffusion coefficient correlations by Wilke and Chang 34 for inorganic, and Young et al. 49 for organic solutes, most aqueous solutions at elevated temperatures attain this coefficient in the range D AB ∈ [1 × 10 −10 , 5 × 10 −8 ] m 2 /s. Consequently, we use this relevant range to investigate the effect of the diffusion coefficient on the solution supersaturation. We aim in this investigation case to examine the qualitative effect of the diffusion coefficient on the supersaturation. To more accurately study a specific solute would require solute-specific models for the diffusion coefficient depending on both solute concentration and temperature.

Case 2: Influence of Laser Pulse Parameters on Supersaturation.
In the second investigation case, we use the 1D framework to examine how the laser pulse energy and the volume of superheated liquid affect the supersaturation of a given solution in a thermocavitation event. For this investigation, we choose the aqueous ammonium sulfate solution, but, in general, the investigation methodology can be applied to study any solution. The results of this investigation case are presented in section 3.4. The laser pulse duration is t las = 6 ns and all material properties, except the solute diffusivity and solubility, are assumed those for pure water. Next, we determine the range of laser pulse energies that is examined in this investigation case and, subsequently, the range of the superheated liquid volume that the laser energy is distributed in.
The laser pulse locally heats the solution due to absorption of the laser irradiation. If the vapor and liquid phases are in thermal equilibrium, the liquid temperature must reach for evaporation to be possible and for the vapor bubble to grow. 50 This condition arises from the Laplace pressure. The vapor bubble studied in this case has an initial radius of R 0 = 0.5 μm, which gives a liquid temperature of at least 416 K for bubble growth to be possible. The liquid temperature is reduced during the laser pulse by thermal diffusion and cooling by evaporation. For a stable growth of the initial bubble, a liquid temperature of 420 K is chosen as the lower temperature limit. This temperature determines the lower limit of the laser pulse energies that we examine in this investigation case. At around 647 K and 22 MPa, water transitions into a supercritical phase and the physical properties change dramatically. This phase is not considered in our numerical framework, and, therefore, the upper limit for the liquid temperature is 610 K, well below the critical temperature. The range (420, 610) K provides the laser energy density limits at the bubble interface, e las (R 0 ), according to We also need the maximum energy density, at the center of the beam, e max , to compute the total laser pulse energy, E las , that produces a desired energy density at the bubble interface. This value is given by eq 21 as Then, the laser pulse energy E las is determined by eq 22 for a given laser beam radius r las . The value of E las is computed so that the energy density at the bubble interface is varied within the range (e R e R ( ) , ( ) las 0 min las 0 max ), for a given laser beam radius. This ensures that also the solution temperature is varied within the range (420, 610) K. The laser beam radius governs the distribution of the laser pulse energy, and the range for this parameter also needs to be defined.
The solution volume that the laser pulse energy is absorbed in can be adjusted by changing the focus of the laser beam. The spatial distribution of the energy should have a three-dimensional character, but, as the vapor bubble grows, the heated liquid region forms a layer adjacent to the bubble interface. This layer becomes thinner as the bubble grows and it is, in the 1D framework, assumed spherically symmetric. Furthermore, we assume a Gaussian irradiation profile of the laser beam given by eq 20, where the value of r las represents the beam radius and determines the radial extent of the heated liquid region.
The focal spot size of the laser beam has a minimum value due to light diffraction that is about 1 μm in experimental setups. 8,11 Therefore, r las,min = 1 μm is the lower limit that we examine in this study.
There is no obvious upper limit for the laser beam radius, but a larger radius requires higher laser energies to produce thermocavitation that may pose a practical limit. To clearly show how the laser beam radius affects the supersaturation at the bubble interface, we choose a maximum laser beam radius 3 orders of magnitude larger than the minimum, r las,max = 1 mm.
By using the 1D numerical framework, we perform simulations across the entire parameter spaces defined in the two investigation Crystal Growth & Design pubs.acs.org/crystal Article cases. The results from these simulations are given in the following section.

RESULTS AND DISCUSSION
Here, we first present the results from our 2D multiphase DNS simulation, described in section 2.1.3, of a laser-induced thermocavitation bubble with experimentally observed crystallization. Afterward, we give the results from the first validation case and the two investigation cases simulated using our 1D framework and described in sections 2.2.1, 2.2.2, and 2.2.3, respectively.
3.1. 2D Simulation Results. A 2D axisymmetric simulation is performed of a laser-induced thermocavitation bubble using a similar setup as in the experiment by Soare et al. 10 We recall that in this experiment crystallization is observed in connection to the thermocavitation event, and that the purpose of the simulation is to test the nucleation hypothesis related to the high solution supersaturation around the bubble during the bubble growth phase.
The solute concentration and the relative supersaturation of the solution at the bubble interface are shown in Figure 3. The saturation mass fraction c sat is computed using the correlation by Daudey 51 (35) that is valid in the range −10 to 400°C. In the experimental study, it is suggested that the nucleation takes place at the maximum rate of evaporation. The evaporation mass fluxes along the r-and z-axes of the computational domain are shown in Figure 4. Here, both curves reach a maximum at around 0.1 μs, and in Figure 3 the relative supersaturation, σ = (c − c sat )/ c sat , reaches above zero at about 0.11 μs. Since nucleation is possible at positive supersaturation values, it seems probable that there exists a close connection between the maximum evaporation rate and the observed crystallization.
The duration of the positive supersaturation is about 1.7 μs, and, after this period, the saturation is reduced to around σ = −0.1 due to diffusion of solute. As the interface liquid continues to cool down, the saturation increases and would reach saturation condition again if the simulation were continued. To assess if the supersaturation peak predicted in the simulation is long enough for primary nucleation to occur, the induction time of the first nucleus is estimated as where J is the nucleation rate per unit volume and V S is the solution volume. 52 By assuming heterogeneous primary nucleation and that the uptake of the growth units by the nucleus is diffusion controlled, the nucleation rate can be estimated from classical nucleation theory according to 23 i k j j j j j y where S = σ + 1 is the supersaturation and A and Δμ are given by  The laser pulse rapidly heats the solution which causes the solubility to increase. This results in an undersaturation that rapidly increases as the solute concentration rises due to solvent evaporation. The supersaturation reaches a maximum of more than 0.13, which is more than can be obtained using evaporative or cooling crystallization under normal conditions. 53,54 Here, D AB = 1 × 10 −8 m 2 /s is the diffusivity of the solute, C 0 = 5 × 10 20 m −3 is the concentration of the foreign particles, k B = 1.381 × 10 −23 J/K is the Boltzmann constant, ν s = 1.35 × 10 −28 m 3 is the partial molecular volume of the solute, ν 0 = 1.24 × 10 −28 m 3 is the molecular volume of the crystalline species, p 0 = 101 300 Pa is the system pressure, and γ = 0.003 J/m 2 is the effective interfacial free energy. 23 A conservative estimate of the thickness of the supersaturated solution around the bubble is ≈1 × 10 −8 m. With this extent, the supersaturated solution volume becomes V S = 1.27 × 10 −19 m 3 . During the supersaturation peak of roughly 1 μs, in Figure  3, the solute concentration, supersaturation, solution temperature, and pressure are approximately c = 0.65, S = 1.1, T = 450 K, and p = 0.9 MPa, respectively. These values into eqs 36−39 give an estimated induction time of t i = 1 × 10 −11 s, which indicates that the duration of high supersaturation is enough for primary nucleation to occur.
Although the primary nucleation is possible, it is not certain if the duration of the supersaturation is enough for the nuclei to grow enough to survive the subsequent period of undersaturation. To predict the correct nucleation and growth kinetics, additional models are needed. These models should be coupled together with the evolution of the solute concentration since the crystal growth collects the available solute. The growth will reduce the solute concentration and the levels of the super-and undersaturation at later stages of the thermocavitation event, which, in turn, affect the following crystal growth or dissolution phases. Therefore, the nucleation and growth kinetics must be predicted accurately. However, as noted in the modeling approach of Soare, 23 the available models for the crystallization kinetics did not provide reliable results at the time scales of the thermocavitation event. In that study, a fitting parameter in the order of 10 5 was needed to match the experimental crystal growth measurements. Therefore, further experiments and development of the nucleation and growth kinetics models, at the relevant time scales, are needed before the models can be used to produce reliable predictions about the crystallization and supersaturation after the primary nucleation has taken place.
The relative supersaturation attains a maximum of more than 0.13, at about 1 μs. This degree of supersaturation is higher than what can be obtained in the aqueous ammonium sulfate solution using evaporative or cooling crystallization under normal conditions. 53,54 If the crystals are nucleated during the bubble growth phase, due to high supersaturation, the simulation clearly indicates that the nucleation occurs in the first few microseconds after the laser pulse. This time frame for nucleation conforms well with the occurrence of the optical disturbances in the experiment which were thought to appear due to nuclei that had grown large enough to become visible.
The liquid temperatures at the bubble interface along the rand z-axes are shown in the top panel of Figure 4. Along the raxis, the superheated liquid region between the bubble interface and the surrounding liquid is cooled by evaporation and thermal diffusion to the surrounding liquid. Here, the thickness of the superheated liquid region is at a minimum, and the cooling rate of the liquid, after the initial fluctuations, is the highest. Along the z-axis, the superheated liquid region is cooled by evaporation and thermal diffusion to the glass slide, which is kept at a constant temperature of 393 K, after the laser pulse. This temperature is 100 K higher than that of the surrounding liquid and results in lower thermal diffusion losses along the z-axis than along the r-axis. The low diffusion heat loss and the thicker region of superheated liquid result in the slowest cooling of interface liquid along the z-axis.
However, the cooling of the interface liquid is mostly due to evaporation, and a heat loss by diffusion is comparably small. Between 0.1 μs and 1 μs, the evaporation rates, shown in the bottom panel of Figure 4, are similar at both axes. Therefore, the cooling rates of the interface liquid at both axes, shown in the top panel of Figure 4, are also similar. The average cooling rate during this period is around 18.9 × 10 6 K/s. The simulation is stopped when most of the superheated liquid is cooled below saturation conditions and the evaporation stops. Since the evaporation governs the increase of solute concentration at the bubble interface, the solution is not expected to reach above saturation conditions again if the simulation is continued.
The evaporation rate at the r-and z-axes represents the two extremes along the bubble interface. These fluxes are similar during the period of high evaporation rate and supersaturation, which indicates that our method, described in section 2.1.1, of taking averaged values across the interface to compute the solute concentration and supersaturation is a reasonable one.
When the evaporation rate decreases and the bubble interface decelerates, the pressure at the z-axis is reduced. This reduction is a consequence of the inertia of the decelerating liquid that flows radially away from the z-axis due to the constraints of the glass slides. The pressure decrease gives a lower saturation temperature at the z-axis than in other parts of the domain. The relatively low liquid cooling rate and the low saturation temperature prolong the evaporation period at the z-axis compared with other parts of the domain. The slow cooling and low pressure are both wall effects that may cause the solution supersaturation to reach higher local values. This indicates that the presence of walls can increase the probability for nucleation and should be taken into account for a precise control of the supersaturation level. With this simulation, however, our aim is to test the plausibility of the nucleation hypothesis, and for that purpose we consider only averaged interface values to obtain a conservative estimate of the supersaturation.
The bubble radius along the r-axis is shown in Figure 5 together with the experimental data. The bubble growth rate is in fair agreement with the experiment, although it is somewhat underestimated. For a comparison during the early bubble growth phase, experimental data sampled with a higher resolution would be needed. The lower growth rate in the simulations can result from an underestimation of the solution temperature after the laser pulse. The 494 K used in this simulation is taken from the experimental study. 23 A higher temperature would cause the evaporation rate to increase, and the bubble to expand faster and give a higher solution supersaturation. This makes the predicted supersaturation maximum of more than 0.13 a low estimate of the experimental case.
To assess if all the correct bubble dynamics is captured, it would be interesting to compare the simulated and measured bubble radius evolutions during the whole bubble lifetime. However, that would require a much longer simulation time than the current 6 μs. Because of the immense computational cost of the 2D simulation, it was not feasible to extend the simulation time to further experimental data points. To reach the present simulation time, the simulation ran on 160 CPUs for two months and required almost 270 000 time steps with a mean time step size of around 2 × 10 −11 s. Nonetheless, as Crystal Growth & Design shown in Figure 3, the period of high supersaturation, that we want to investigate in this work, is fully captured in the current simulation time. The high computational cost is indeed the main argument for using the 1D model to perform the more extensive parameter investigation cases.
In summary, our 2D simulation results agree with the nucleation hypothesis which we considered plausible for the studied thermocavitation event.
3.2. Results from the First Validation Case (1D Model). In our first validation case of the 1D model, we compare the predicted bubble radius evolution with an experimentally observed laser-induced thermocavitation event reported by QuintoSu et al. 48 The comparison is shown in Figure 6, and we note that a good agreement between the two is found, which indicates that the 1D framework predicts realistic evaporation rates. This validation case is described in section 2.2.1. The results from the second validation case of the 1D model can be found in the Supporting Information 2.2.
In addition to the bubble evolution, we present some other interesting fluid dynamic aspects of the thermocavitation event, as predicted by the 1D framework. Figures 7, 8, and 9 show temperature, pressure, and velocity profiles, respectively, at three instants of the simulation.
The first instant is taken at 0.5 ns, shortly after the laser pulse. The interface liquid is superheated and far from saturation conditions, which results in an intense evaporation. The saturation pressure of the superheated liquid is around 10 MPa, and the vapor pressure is rapidly approaching that value. The increased vapor pressure produces a vapor temperature of several thousand Kelvin. Also, high velocity waves propagate in the bubble due to the evaporation and pressure fluctuations.  Soare. 23 The temporal resolution of the experimental data is 4 μs and can therefore not be used for a direct comparison of the early bubble growth phase. The simulation underpredicts the growth rate, which can indicate that the solution temperature after the laser pulse should be higher than the 494 K estimated in the experimental study. A higher temperature would result in a higher solvent evaporation rate, a faster bubble expansion, and a higher interface concentration. Figure 6. Bubble radius evolution from our 1D framework compared with experimental data. 48 The estimated uncertainties of the measurements are ±0.5 μs and ±2 μm. In the experiment, the laser-induced thermocavitation bubble was formed with a glass wall on one side and unbounded elsewhere, while the simulation uses an entirely unbounded domain. The calculated and measured bubble radius evolutions are in a reasonable agreement, which ensures that realistic bubble dynamics and evaporation rates are predicted by the 1D framework. Figure 7. Temperature profiles at three instants in the 1D simulation of an experimentally observed thermocavitation bubble. 48 The bubble has an initial radius of 0.5 μm, and the laser beam radius is 11 μm. The black lines represent the position of the bubble interface, and the arrows indicate the interface velocity direction. Crystal Growth & Design pubs.acs.org/crystal Article At 0.1 μs, the vapor pressure is almost constant and reduced to the lower saturation pressure of the liquid that is cooled by evaporation and the surrounding liquid. The vapor temperature also decreases due to cooling by the liquid and the pressure reduction. Still, the vapor pressure is above the surrounding atmospheric pressure, which results in the bubble growth. At this instant, the liquid velocity at the interface reaches more than 40 m/s, while the velocity of the vapor at the interface is around 10 m/s due to the evaporation and the density difference of the phases.
After about 1.3 μs, the evaporation stops since the interface liquid is sufficiently cooled by the latent heat of evaporation and thermal diffusion. The bubble continues to expand for a few microseconds due to liquid inertia, which reduces the bubble pressure below saturation conditions and results in vapor condensation. These phenomena produce a minimum vapor pressure of around 5 kPa. This pressure is much lower than the surrounding atmospheric pressure which causes the bubble to decelerate and then shrink.
At the last instant, at 11 μs, the bubble is close to collapse, and the interface approaches sonic velocities. Since we do not include compressibility effects in the liquid phase, all aspects of the collapse are not captured, but the results clearly show the rapid increase of pressure and temperature of the vapor as the inertia of the high velocity liquid compresses the bubble. Comprehensive simulations of the collapse phase in other studies demonstrated outgoing liquid shock waves, with temperatures reaching almost 10 000 K and pressures greater than 7 GPa. 44 3.3. Results from Case 1 (1D Model): Influence of Solution Properties on Supersaturation. In this section, we present the results obtained by our 1D framework from our first investigation case of the nucleation hypothesis. Here, we study the effects of the solubility and the solute diffusivity on the interface solution supersaturation. The simulations are performed with the parameters (C 1 , n, D AB ) varied within the specified ranges that correspond to realistic and industrially relevant conditions. From each simulation, we determine the maximum degree of interface supersaturation during the entire thermocavitation simulation. The results from these simulations are shown in Figure 10.
Three values of the exponent n, in the temperaturedependent solubility curve, eq 24, are used. At a given temperature, a higher value of n results in a higher solubility and hence a lower supersaturation. The solution temperature and bubble evolution are the same for all simulations in this investigation case since we assume pure water material properties and use a fixed set of laser pulse parameters. Therefore, a higher value of n results in lower supersaturation levels.
The range of the solubility parameter C 1 is chosen so that the supersaturation σ does not change significantly at higher or lower values of C 1 , as discussed in section 2.2.2. For solutes that are highly soluble at 0°C and have a flat solubility curve, the solubility parameter C 1 is small. Conversely, less soluble solutes with steep solubility curves have a higher value of C 1 . For a solution heated by a laser pulse, the lower the parameter C 1 , the lower the increase of the solubility relative to the solute concentration. The lower relative increase of the solubility results in a higher supersaturation σ. For solutions with higher values of C 1 , the increase of the solubility, relative to the solute  Three values of n are used that represent different types of behavior of the temperature-dependent solubility curve. Higher values of n give a higher solubility at a given temperature and therefore lower supersaturation. For visualization purposes, the supersaturation is limited to the range σ ∈ [0.01, 10]. The lowest value of all simulations is σ = 0, which is the initial condition, and the highest value is σ = 35, obtained with the lowest values of D AB and C 1 . The symbol "×" represents the aqueous ammonium sulfate solution with a maximum supersaturation of more than 0.9. The symbol "+" represents an aqueous potassium chlorate solution that does not reach the above explained saturation conditions during the entire simulation.
Crystal Growth & Design pubs.acs.org/crystal Article concentration, is higher, which results in a lower degree of supersaturation.
The evaporation of the solvent increases the solute concentration at the interface, while the resulting concentration gradient induces solute diffusion away from the interface. The net effect of those two phenomena determines the concentration of solute at the bubble interface. Since the evaporation rate is the same for all the simulations in this investigation case, we can study the influence of the diffusion mass flux independently. A lower value of the solute diffusivity, D AB , gives a lower diffusion mass flux, for a given concentration gradient. The lower diffusion mass flux results in a higher interface solute concentration and thereby a higher degree of supersaturation. Higher values of D AB give a higher diffusion mass flux, lower concentrations of solute, and lower supersaturation at the bubble interface.
The diffusivity value corresponding to the aqueous ammonium sulfate solution is indicated in Figure 10 by the symbol "×". The ammonium sulfate is highly soluble, and the solubility curve is approximately linear with a low slope that results in a low solubility parameter of C 1 ≈ 0.044. In the thermocavitation event studied in this section, the ammonium sulfate solution attains a maximum interface supersaturation of more than 0.9. This degree of supersaturation is, like in the 2D simulation of section 3.1, higher than what can be obtained using evaporative or cooling crystallization under normal conditions. Also, the estimated value for an aqueous potassium chlorate solution is shown in Figure 10 by the symbol "+". This salt solution has a significantly higher solubility parameter of C 1 ≈ 1.98, which indicates a steep solubility curve. After the heating by the laser pulse, the solubility is higher than the solute concentration reaches during the entire thermocavitation event.
The simulation results clearly indicate that nucleation resulting from high interface supersaturation is highly dependent on the solution properties and not trivial to obtain for all solutions. If the solubility can be approximated by eq 24, the produced maps in Figure 10 can be used as guidelines to determine whether a specific solution may reach conditions favorable for primary nucleation given the set of realistic laser pulse parameters. In general, the methodology used to produce the maps in Figure 10 can be modified to investigate any type of solution.
3.4. Results from Case 2 (1D Model): Influence of Laser Pulse Parameters on Supersaturation. Here, we present the results from our second investigation case where we study the effects of the laser pulse energy and the superheated liquid volume on the interface solution supersaturation. In this case, the absorbed laser pulse energies are chosen so that the energy density in the interface liquid e las (R 0 ) produces liquid temperatures within the range (420, 610) K. The reader is referred to section 2.2.3 for a detailed description of how the laser pulse energies are determined. The laser beam radius r las determines the spatial distribution of the laser energy density, according to eq 21. In this investigation case, we have chosen laser beam radii in the range (1, 1000) μm.
The maximum relative supersaturation obtained from the simulations is shown in Figure 11. In these simulations, we use the aqueous ammonium sulfate solution. Initially, the solute concentration is at saturation conditions of around 43 wt% at 20°C. When the superheated water is evaporated, the solute concentration increases, and for high laser energy densities the evaporation is so intense that the 1D framework predicts concentrations of more than 100 wt%. These results all lie above the dashed line in Figure 11 and are possible in the numerical framework since we treat the solute concentration as a passive scalar, assume material properties of pure water and do not include the effects of crystallization. The formation of a solid crystal phase would reduce the solute concentration of the solution, thus preventing the concentration to reach such high values. In addition, the water activity should decrease as the solute concentration increases. A lower water activity reduces the evaporation rate and thereby also limits the increase of the solute concentration. The effect of the water activity is, however, not included since we assume the material properties of pure water. By not accounting for these limiting effects on the concentration, and by treating the concentration as a passive scalar, it is possible for the concentration to reach nonphysical values at the high energy densities. To more accurately study a specific solution would require material property models that also depend on the concentration. Still, we argue that the results below the dashed line are the physically realistic estimates of the supersaturation and that they clearly show how the investigated laser parameters affect the supersaturation of the solution at the bubble interface.
A higher laser energy density e las (R 0 ) produces a higher liquid temperature and, as a consequence, faster evaporation of the solvent. When the solvent is evaporated more rapidly, the concentration of solute reaches higher values, and the supersaturation of the solution increases. Similarly, for a given energy density, a larger laser beam radius also produces a higher solution supersaturation. The larger radius superheats a greater volume of liquid, and, although the liquid temperature governs the evaporation rate, the volume of the superheated liquid determines the duration of the evaporation. A longer period of solvent evaporation gives a higher concentration of solute and thereby also a higher solution supersaturation. Figure 11. Maximum interface supersaturation of the aqueous ammonium sulfate solution during the 1D thermocavitation bubble simulations of our second investigation case. In this case, we study the influence on the solution supersaturation by the laser beam radius, r las , and the laser energy density at the bubble interface, e las (R 0 ). For visualization purposes, the supersaturation is limited to the range σ ∈ [0.01, 1]. The simulation cases above the dashed line attain interface solute concentrations of more then 100 wt%, which is numerically possible but not physical. Therefore, the results from these simulation cases are not included. Below the dashed line, the results are physical estimates of the solution supersaturation. The bullet symbol (•) represents the laser pulse parameters equivalent to those used in the 2D simulation of the same ammonium sulfate solution. The results shown in Figure 11 clearly demonstrate that greater laser energy densities and laser beam radii produce higher solution supersaturation levels. This tendency is the same for all solutions, although the laser parameters, which produce a certain degree of supersaturation, vary depending on the solution properties. Simulations with other solubility parameters C 1 are performed and indeed show the same trends. In general, the methodology used in this investigation case can be adopted for any solution to investigate possible degrees of supersaturation and to obtain laser pulse parameters that produce conditions favorable for primary nucleation.
Material properties different than those of pure water would affect the degree of supersaturation. Specifically, the effect of the water activity would lower the evaporation rate and thereby reduce the increase of solute concentration and supersaturation. We expect, however, that the general findings and trends of our investigations hold and that the results would not change significantly for most other aqueous solutions.
For the sake of a comparison of the predicted supersaturation levels obtained by the two numerical frameworks used in this work, the black bullet symbol (•) in Figure 11 represents equivalent laser pulse parameters as those used in the 2D simulation case. In the 2D simulation, the ammonium sulfate solution is heated uniformly within a cylindrical region, while in the 1D simulations the solution is heated using a spherical Gaussian power distribution. Strictly speaking, these differences do not allow for a direct comparison of the results. Yet, the equivalent laser pulse parameters lie in a region of the parameter space with supersaturation values close to the 2D simulation result of σ = 0.13. The comparison shows that the 1D and 2D frameworks predict similar supersaturation levels at equivalent laser parameters despite different laser energy distributions. This indicates that our assumption of a spherically symmetric energy distribution in the 1D framework is reasonable. A more detailed discussion of the predicted interface supersaturation by the two numerical frameworks is presented in the Supporting Information 2.3, and we note that a good agreement between the two is found.
The laser energy density e las (R 0 ) determines the solution temperature at the bubble interface after the laser pulse. This is the maximum solution temperature that is obtained in the liquid during the thermocavitation event except for possible peaks at the interface during the bubble collapse. For solutes that are sensitive to thermal degradation it is therefore favorable to suppress the value of e las (R 0 ) that is necessary for nucleation of a given solution. Interestingly, the results of Figure 11 show that high supersaturation values can be obtained at relatively low values of e las (R 0 ) if a sufficient laser beam radius r las is used. A larger laser beam radius produces a greater volume of superheated liquid and a longer period of evaporation. The prolonged evaporation period gives higher solute concentrations at the interface without the need for higher evaporation rates by increased solution temperatures. Thus, to achieve high supersaturation at a relatively low temperature, the superheated liquid volume should be large and the laser energy density minimized.

CONCLUSIONS
We tested and confirmed the plausibility of the hypothesis that high supersaturation is produced, due to solvent evaporation, during the early growth period of laser-induced thermocavitation bubbles. Using a 2D-axisymmetric multiphase DNS simulation of a thermocavitation event combined with an experimentally observed crystallization, we obtained a conservative estimate of the relative supersaturation reaching above 0.13. This degree of supersaturation is higher than what can be obtained using conventional crystallization techniques under normal conditions and indicates that the nucleation hypothesis is viable.
In an effort to deal with a prohibitive cost of 2D DNS simulations and, consequently, the lack of possibility to provide a suitable parameter space for the analysis, we developed a 1D numerical framework that is able to accurately reproduce the entire thermocavitation event, from the laser pulse to the bubble collapse phase. Using this framework, we conducted two extensive investigation cases that examined how the solute solubility, solute diffusivity, laser pulse energy, and volume of superheated liquid affect the supersaturation of the solution around the thermocavitation bubble.
In the first investigation case, we have shown that the degree of supersaturation is highly dependent on both the solubility and the diffusivity of the solute. For the same laser pulse parameters, an ammonium sulfate solution attained a maximum supersaturation of more than 0.9, while a potassium chlorate solution did not reach above saturation conditions during the entire thermocavitation event. The simulation results showed that high interface supersaturation was not readily obtained for all types of solutions in a typical thermocavitation event. The produced maps, in section 3.3, can be used as a guideline to identify if a specific solution is probable to reach conditions favorable for crystallization, given the set of realistic laser pulse parameters used in the simulations.
We showed in the second investigation case that higher values of the laser pulse energy and the superheated liquid volume resulted in a higher solution supersaturation. However, not all combinations of the examined laser pulse parameters produced supersaturation of the studied ammonium sulfate solution. A careful selection of the laser pulse parameters is necessary to obtain the desired degree of supersaturation. Again, the produced map, in section 3.4, may be used as a guideline to find an appropriate set of laser pulse parameters. In addition, we showed that the combination of a large volume of superheated liquid and a low laser energy density can be used to obtain high supersaturation at a relatively low temperature. This type of laser parameter is favorable for solutes sensitive to thermal degradation.
In conclusion, all simulation results presented in this study clearly show that the laser-induced thermocavitation method can produce high solution supersaturation under a range of realistic and industrially relevant conditions. The numerical frameworks and insights provided in this work can be used to investigate and design appropriate laser-induced nucleation setups for obtaining high solution supersaturation and conditions favorable for crystallization.
The proposed nucleation hypothesis and our conclusions could be further supported by combining the type of simulations performed in this work with new experiments, as outlined in section 1. Such investigations may elucidate connections between the degree and duration of solution supersaturation and the crystallization.