Simulations of shock waves and cavitation bubbles produced in water by picosecond and nanosecond laser pulses

We compare numerical simulations of bubble dynamics in water with experiments performed at the Medizineisches Laserzentrum Lubeck. Spatial and temporal features of the laser beam were modeled. Plasma growth was predicted using a moving breakdown model. We compare the measured and calculated positions of the shock front and the bubble wall as a function of time after optical breakdown in water. Nd:YAG laser pulses of 30-ps 1-mJ and 6-ns 10-mJ were simulated. We have extended previous work in which picosecond deposition was modeled at temporally instantaneous and spatially uniform.


INTRODUCTION
Bubble dynamics and its associated shock waves, acoustic waves, and jetting, are important mechanisms in laser medical applications. The dynamics produce both intended medical effects and undesirable collateral effects. Multidimensional Eulenan compressible hydrodynamics codes are capable of modeling the shock waves and highly distorted flows associated with these applications. Combined with appropriate laser absorption models, material equations of state (EOS), tissue dynamics and tissue failure models, such codes can aid in understanding and optimizing clinical protocols in which bubble dynamics are an important consideration. Vogel et al.24 used time-resolved photographic methods to investigate shock wave and cavitation bubble expansion after optical breakdown in water using Nd:YAG laser pulses of 30-ps and 6-ns duration (full width at half-maximum) with energies typical of medical applications. They measured the position of the shock front and bubble wall position as a function of time with temporal resolution '4-ns and spatial resolution -.-4-m. These high-resolution measurements provide excellent data for quantitative comparison with numerical simulations. Selected photos from the Medizinisches Laserzentrum LUbeck (MLL) experiments are reproduced in Figs. 1 and 2. The laser light is incident from the right.
In previous work by Chapyak et al. the Los Alamos National Laboratory code MESA-2D was used to examine the MLL picosecond experiments assuming that the energy deposition was temporally instantaneous and spatially uniform. These assumptions worked well because there is little dynamic motion of the target material during energy delivery for the short pulse width. The excellent agreement of the calculations with the experimental data serves to validate both the numerical hydrodynamics of the MESA-2D code and the water equations of state used in the simulations. The principal purpose of the present calculations is to extend the previous work to include the 6ns MILL experiments. For these experiments there is significant dynamic response of the target material during the laser pulse and we could no longer make the assumption that the energy deposition is temporally instantaneous and spatially uniform. In order to accomplish the simulation we modeled the temporal and spatial features of the laser beam, plasma ignition. plasma growth, and the deposition of the laser energy in the plasma. Relatively straightforward models were employed which are easy to implement and use but which ignore the detailed plasma physics. Plasma ignition and growth is predicted using a moving breakdown model while the laser energy deposition in the plasma is treated as an exponential function.

MESA-2D ENERGY DEPOSITION MODEL
In order to treat the temporal and spatial aspects ofthe laser energy deposition it was necessary to add a new model to the MESA-2D user-defined subroutines. The model accounts for the pertinent aspects of the laser experiments presented in Ref. 2; i.e., the beam focusing angle, Gaussian radial energy distribution and the experimentally determined energy rate.
The laser beam energy is tracked using constant-energy beam segments based on the Eulerian mesh at the beam focal plane and includes the Gaussian radial energy distribution. These "constant energy" segments are projected back toward the laser and used to follow the incoming energy. The model follows the beam out to 1.5 times the lie2 radius and thus tracks 98.9 % of the beam energy under the assumption that it is Gaussian.
Ionization is predicted using the moving breakdown model discussed in Refs. 3 and 4. It assumes that breakdown occurs independently at each location as soon as the irradiance of the electrical field surpasses the threshold value. At each time step in the calculation the energy irradiance is calculated as a function of time, and spatial location in the beam, with the result checked against the ignition threshold.
Energy deposition in the plasma is treated as exponential, as a function of distance from the leading edge of the plasma. No energy is deposited in cells that are below the breakdown threshold and laser energy is either absorbed or lost via transmission. Only minimal energy reflection was measured in the experiments4 and thus none is included in the model. The nodel keeps track of the energy absorbed and lost.
References 3 and 4 report the ignition thresholds and absorption coefficients that were derived from experimental data. These were used as input parameters for the energy deposition model and are listed in Table I. Both values are treated as constants. The temporal evolution of the laser power PL is described as a sin2 function with duration r (full-width at half-maximum), and total duration 2 v pp0sjfl2(tJ Table I. Laser parameters derived from experiments 30-ps laser 6-ns laser Breakdown irradiance threshold 6.3x10' W/cm2 5. lxlO'° W/cm2 Exponential absorption 360 cm' 900 cm' Several alternative energy deposition options are included in the general model. All use the moving breakdown model to predict ionization but differ in the way in which the energy is distributed within the plasma. Options include uniform energy deposition in the plasma as a function of cell volume, or mass. Another option deposits energy as a function of volume in each beam segment but retains the Gaussian radial energy distribution across the beam segments. These options allow investigation, in a general sense, of the affects of energy redistribution within the plasma.

CONDITIONS FOR THE SIMULATIONS
The pertinent laser parameters for the experiments, as defined in Ref. 2, are listed in Table II. We used the same LANL SESAME5'6 tabular water equation of state (EOS) that was used in the previous work. Identified as SESAME #7150, it is an equilibrium, Maxwell construction EOS that includes the effects of both phase changes and ionization. The simulations do not include the effects of viscosity or surface tension.

COMPUTATIONAL RESULTS AND COMPARISONS TO EXPERIMENT
As discussed above, the primaiy purpose of the present work was to simulate the 6-ns experiment where there is significant material motion during the laser pulse. In addition to the long pulse-width simulation we reran the 30-ps case using the new model to include the time dependent plasma dynamics and energy deposition.
The computed and measured maximum radial shock wave and bubble wall locations for the 6-ns, lO-nU case are plotted as a function of time in Fig. 3. The shock wave location for the simulation was defined as the location of the 0. 1 kbar pressure contour while the bubble wall was defined as coincident with the 0.98 g/cm2 density contour. Agreement is reasonably good with the velocity of both the shock and bubble wall being a little higher than the experiment at the end of the calculation. The discrepancy of the starting position could be caused by differences between what is apparent in the photos versus what was selected as the edit in the code simulation.
Selecting the element size for the computational mesh is always a trade off, particularly for multidimensional simulations. Finer zones result in finer resolution at the cost of computer time that can result in run times that are cumbersome (long turn around) or prohibitive (cost or resources). The present calculation is a good demonstration of this and of how the velocity of the shock can be adversely affected by the mesh. As discussed previously, in order to minimize the total number of elements in the problem the cell size outside of the plasma volume is increased using a multiplicative factor. The original factor for this problem was 1. 1 . The resulting mesh tended to spread out the wave front resulting in apparent shock velocities that were unrealistically high. Switching to a rezone ratio of 1 .03 more than doubled the computational time but resulted in more realistic shock velocities. The logical next step of further refining the mesh has not been done.  The energy distribution computed by MESA-2D is plotted in Fig. 4. The initial energy in the simulation is defined as zero for the purpose of the plot, it is not zero in the calculation. Figures 5and 6 are the computed pressure and density as a function of the radius, at the plane of the maximum shock wave and bubble wall radius. The pressure profile at 90-ns shows considerable spreading of the shock front which is probably a function of the increasing mesh size even with the revised rezone ratio.
Hydrocodes normally must spread the shock front over several cells to maintain calculational stability. The calculations can not handle an infinitely steep shock front. The density plot shows a high density region following the shock wave with a sharp decrease in density at the bubble wall.
Figures 7 and 8 are contour plots of the bubble wall (density = 0.98 g/cm3) and the shock front (jressure = 0. 1 kbar). These are roughly the computational equivalent of the experimental photos in Fig.  1 . There are similarities but also a number of obvious differences in shape. The steps in the bubble contour in Fig. 7 are related to the beam segments used in the energy deposition calculation. The numerical instability or noise in the shock contour in Fig. 8 is essentially inherent in the computational process. It is not a function of the energy deposition model and does not affect the results. It is possible that refining the mesh would eliminate it but since the results are not affected, and considering the increased run time, this was not pursued.
A limited number of calculations were run to examine the sensitivity of the results to various input parameters and the optional energy deposition models. None of calculations resulted in a significant improvement however there are still a number of options to be investigated.  The computed and measured maximum radial shock wave and bubble wall location for the 30-ps, 1-mJ case are plotted as a function of time in Fig. 9. Again the shock location for the simulation was defined as the location of the 0. 1 kbar shock front while the bubble wall was defined as being coincident with the 0.98 g/cm2 density contour. Agreement is quite good for the position and velocity of the shock front. However, there is rather poor agreement between measured and computed bubble wall position and velocity.
The reason for this discrepancy remains unclear. The model works for the 6-ns case and checks of the model as incorporated in the code have not revealed any problems. A limited input parameter study, including the various energy deposition options, did not result in marked improvement. The simulation plotted in Fig. 9 used the uniform-mass energy deposition option that distributes energy uniformly throughout the plasma volume as a function of cell mass at the time the energy is deposited. It provided marginally better results than the standard exponential deposition option. The lack of agreement may indicate that our relatively simple model is not an adequate representation of the laser-matter interactions in this experiment.
The energy distribution computed by MESA-2D is plotted as a function of time in Fig. 10. Figures 1 1 and 12 present the computed pressure and density as a function of the radius at the plane of the maximum shock wave and bubble wall radius. Again the pressure profile at 90-ns shows considerable spreading of the shock front. The density plot shows the low density region starting at the outside of the plasma region and working its way inward as the pressure is relieved. Figures 13 and 14 are contour plots of the bubble wall and the shock front. These are roughly the computational equivalent of the experimental photos in Fig. 2. In this case the differences between the computed contours and the photos are more pronounced. Most obvious is the lack of a down-stream tail in the calculated contours. The tail is probably due to self-focusing in the medium in conjunction with self-defocusing inside the plasma. This cannot be modeled with the present energy deposition scheme. In the calculation, essentially all of the laser energy is absorbed in the plasma from the focal plane forward with none left to cause ionization down stream from the focal plane. The steps in the bubble contour in

CONCLUSIONS
We have simulated the dynamics produced by both 6-ns 1O-mJ and 30-ps 1-mJ laser pulse plasma breakdown in water. The computed shock wave and bubble wall locations for the 6-ns case are in reasonable agreement with the measurements made at MLL. For the 30-ps case the shock wave predictions agree quite well with the measurement but there is a problem with the bubble size and bubble wall velocity. It is not clear if the discrepancy can be corrected by further work with the existing model or