1. Introduction
Laser-induced cavitation plays an important role in many fields, especially in laser materials processing in liquid environments (Barcikowski et al. Reference Barcikowski, Plech, Suslick and Vogel2019; Kanitz et al. Reference Kanitz, Kalus, Gurevich, Ostendorf, Barcikowski and Amans2019) and in biomedicine and biophotonics, where it enables precise surgery on cells and within transparent tissues (Vogel et al. Reference Vogel, Schweiger, Frieser, Asiyo and Birngruber1990; Juhasz et al. Reference Juhasz, Loesel, Kurtz, Horvath, Bille and Mourou1999; Koenig et al. Reference Koenig, Riemann, Fischer and Halbhuber1999; Tirlapur & Koenig Reference Tirlapur and Koenig2002; Yanik et al. Reference Yanik, Cinar, Cinar, Chisholm, Jin and Ben-Yakar2004; Vogel et al. Reference Vogel, Noack, Huettman and Paltauf2005; Chung & Mazur Reference Chung and Mazur2009; Palanker et al. Reference Palanker2010; Stevenson et al. Reference Stevenson, Gunn-Moore, Campbell and Dholakia2010; Hoy et al. Reference Hoy, Ferhanoglu, Yildirim, Kim, Karajanagi, Chan, Kobler, Zeitels and Ben-Yakar2014). It involves localized energy deposition by laser pulses that results in an explosive phase transition of the target material causing the bubble expansion. That is accompanied by the emission of an acoustic transient, which often evolves into a shock wave. The bubble expansion opens a cavity in the liquid medium, and its collapse is again accompanied by acoustic transient emission. Shock waves and cavitation are inevitably linked to laser ablation in liquids and to laser surgery. They may contribute to the desired effect but can also evoke undesired side effects, especially when very gentle cuts or perforations in biological environments are desired. Therefore, it is of great interest to systematically explore the dependence of laser-induced cavitation phenomena on laser parameters and bubble size. Moreover, the topic is of general physical interest as laser-induced cavitation involves nonlinear optical, acoustical and hydrodynamic phenomena producing extreme states of matter.
Localized deposition of laser light energy into the bulk of water or transparent biological media relies on plasma formation by optical breakdown of the liquid and is associated with high volumetric energy densities, temperatures and pressures (Vogel et al. Reference Vogel, Nahen, Theisen and Noack1996b; Lauterborn & Vogel Reference Lauterborn and Vogel2013). Therefore, bubbles often exhibit large-amplitude oscillations. With tight focusing and an isotropic environment, almost perfectly spherical bubbles can be produced that match the conditions assumed in theoretical models of spherical bubble dynamics.
In the past, experimental studies of laser-induced cavitation in water have mostly been performed on millimetre-sized bubbles, where time scales are long enough to enable recording of the dynamics via high-speed photography (Benjamin & Ellis Reference Benjamin and Ellis1966; Lauterborn & Bolle Reference Lauterborn and Bolle1975; Tomita & Shima Reference Tomita and Shima1986; Vogel, Lauterborn & Timm Reference Vogel, Lauterborn and Timm1989; Vogel, Busch & Parlitz Reference Vogel, Busch and Parlitz1996a; Philipp & Lauterborn Reference Philipp and Lauterborn1998; Baghdassarian, Tabbert & Williams Reference Baghdassarian, Tabbert and Williams1999; Brujan et al. Reference Brujan, Nahen, Schmidt and Vogel2001; Lindau & Lauterborn Reference Lindau and Lauterborn2003; Brujan & Vogel Reference Brujan and Vogel2006; Obreschkow et al. Reference Obreschkow, Tinguely, Dorsaz, Kobel, de Bosset and Farhat2011; Obreschkow et al. Reference Obreschkow, Tinguely, Dorsaz, Kobel, de Bosset and Farhat2013; Reuter et al. Reference Reuter, Gonzalez-Avila, Mettin and Ohl2017; Sagar & el Moctar Reference Sagar and el Moctar2020; Podbevsek et al. Reference Podbevsek, Lokar, Podobnikar, Petkovsek and Dular2021). Here, surface tension and viscosity play only a minor role and the bubble dynamics exhibits self-similar features over a large range of bubble sizes (Plesset & Prosperetti Reference Plesset and Prosperetti1977; Lauterborn & Kurz Reference Lauterborn and Kurz2010). Laser surgery of biological tissues and, particularly, cells goes along with much smaller bubble sizes in the micro- and nanometre range and shorter oscillation times (120 ns for a bubble with 1 μm maximum radius), which requires faster experimental techniques than high-speed photography. The dynamics of such small bubbles is influenced by surface tension and viscosity because their contributions to the bubble wall pressure scale inversely with the bubble radius R (Lauterborn & Kurz Reference Lauterborn and Kurz2010). Moreover, the laser pulse duration can become a significant part of the bubble oscillation time, especially when nanosecond (ns) laser pulses are used and the bubble expansion starts already during energy deposition. This establishes a need for a systematic investigation of the changes in bubble dynamics with decreasing bubble size, both in biological media and in water. In doing so, it is of the utmost importance to establish an energy balance tracing the partitioning of absorbed laser energy into vaporization, shock wave emission, bubble formation, viscous damping and condensation. Changes in the influence of viscosity and surface tension with decreasing bubble size will strongly affect the partitioning. An energy balance will thus help us to understand the changing dynamics for the transition from micro- to nanocavitation, and it will elucidate the mechanisms governing cell and tissue surgery, and the accompanying side effects.
In this paper, we present a set of experimental and simulation tools that enables a systematic study of the parameter dependence of bubble dynamics including the determination of peak pressures upon bubble generation and collapse, tracking of bubble oscillations and shock wave emission and a complete energy balance. Simulations are based on the well-established Gilmore model (Gilmore Reference Gilmore1952) that is extended by a term covering the initial shock-driven acceleration of the bubble wall in second-order approximation and by tracking the partitioning of absorbed laser energy into vaporization, bubble and shock wave energy and dissipation through viscosity and condensation. Our simulations rely on very few parameters that need to be determined experimentally. These are the size of the laser-produced plasma, which is translated into the start radius R 0 for the simulations, and the duration of the subsequent bubble oscillations, Tosci, where i denotes the number of the oscillation. Numerous studies have confirmed that modelling predictions on R(t) match experimental data for the initial bubble oscillations very well (Lauterborn Reference Lauterborn1974; Müller et al. Reference Müller, Bachmann, Kroninger, Kurz and Helluy2009; Kroeninger et al. Reference Kroeninger, Koehler, Kurz and Lauterborn2010; Obreschkow et al. Reference Obreschkow, Tinguely, Dorsaz, Kobel, de Bosset and Farhat2013). Therefore, the bubble dynamics can be characterized by numerical simulations if the above mentioned benchmark data on plasma size and bubble oscillation times are available.
The plasma size can be determined by time-integrated photography of the plasma luminescence or by time-resolved photography of the optical breakdown region, if the luminescence is too weak (Vogel et al. Reference Vogel, Busch and Parlitz1996a; Schaffer et al. Reference Schaffer, Nishimura, Glezer, Kim and Mazur2002; Venugopalan et al. Reference Venugopalan, Guerra, Nahen and Vogel2002). Bubble oscillation times can be determined with high temporal resolution through single-shot measurements detecting the forward scattering signal from a continuous wave (cw) probe laser beam (Vogel et al. Reference Vogel, Linz, Freidank and Paltauf2008). Various probe beam detection schemes with different beam diameters at the bubble position and different angles of scattered light collection have been employed for bubble monitoring (Barber et al. Reference Barber, Hiller, Lofstedt, Putterman and Weninger1997; Matula Reference Matula1999; Gompf & Pecha Reference Gompf and Pecha2000; Weninger, Evans & Putterman Reference Weninger, Evans and Putterman2000; Schaffer et al. Reference Schaffer, Nishimura, Glezer, Kim and Mazur2002). In the present paper, we use the confocal scheme introduced by Vogel et al. (Reference Vogel, Linz, Freidank and Paltauf2008) in which the probe beam is collinear with the pump beam producing the bubble and the focus locations of both beams coincide. When the DC background is removed by AC coupling of the photodetector, very small oscillations <1 nm in the late phase of the bubble lifetime can be detected, which is not possible photographically. This way, we were able to precisely determine bubble oscillation times through single-shot measurements, and we could trace the transition from the nonlinear large-amplitude oscillations immediately after laser-induced bubble generation to the linear small-amplitude oscillations of the long-lived residual gas bubble.
The experimental generation of highly spherical bubbles requires tight focusing of the pump laser pulses in order to guarantee the formation of compact plasma driving the bubble expansion (Venugopalan et al. Reference Venugopalan, Guerra, Nahen and Vogel2002; Vogel et al. Reference Vogel, Linz, Freidank and Paltauf2008; Obreschkow et al. Reference Obreschkow, Tinguely, Dorsaz, Kobel, de Bosset and Farhat2013; Sinibaldi et al. Reference Sinibaldi, Occhicone, Alves Pereira, Caprini, Marino, Michelotti and Casciola2019). Additionally, buoyancy effects that could lead to a movement out of the probe beam focus and to jet formation upon collapse (Benjamin & Ellis Reference Benjamin and Ellis1966; Zhang et al. Reference Zhang, Cui, Cui and Wang2015) must be avoided. Previous high-speed photographic studies on millimetre-sized bubbles eliminated buoyancy by investigating the bubble dynamics in a falling apparatus (Benjamin & Ellis Reference Benjamin and Ellis1966; Blake & Gibson Reference Blake and Gibson1987) or under zero gravity conditions during parabolic flights (Obreschkow et al. Reference Obreschkow, Tinguely, Dorsaz, Kobel, de Bosset and Farhat2011, Reference Obreschkow, Tinguely, Dorsaz, Kobel, de Bosset and Farhat2013). We limit buoyancy effects by restricting investigations to small bubbles with a maximum radius below 100 μm that are produced by focusing the pump laser pulse through a long-distance water-immersion microscope objective with large numerical aperture (NA). Tight focusing ensures the formation of spherical bubbles by the pump laser beam and provides a high sensitivity of the detection of bubble oscillations down to small bubble sizes.
For the modelling of laser-induced spherical bubble dynamics and shock wave emission, liquid compressibility must be considered. The model by Keller & Miksis (Reference Keller and Miksis1980) does this assuming a constant sound velocity in the liquid, whereas the Gilmore model uses the local pressure conditions at the bubble surface. As a consequence, the Keller—Miksis equation suffers from the unphysical behaviour that acceleration and pressure differences have opposite signs when the bubble wall velocity exceeds the sound velocity in the liquid (Prosperetti & Hao Reference Prosperetti and Hao1999). Gilmore's approach avoids this problem and enables us to follow shock wave propagation into the surrounding liquid even under extreme conditions of optical breakdown, where bubble wall velocities well above 1500 m s−1 are observed (Vogel et al. Reference Vogel, Busch and Parlitz1996a). We complemented it by an automatized procedure for determining the shock front position that facilitates tracking of the pressure decay and energy dissipation at the shock front.
Most previous models of bubble dynamics and shock wave emission started with the expanded bubble and focused on its collapse and rebound, whereas the present model includes also the initial expansion phase. It covers the continuous increase of the driving force for bubble expansion during the laser pulse (following Vogel et al. Reference Vogel, Nahen, Theisen and Noack1996b) and the contribution of the particle velocity behind the detaching shock front producing a jump start of the bubble wall velocity. Gilmore (Reference Gilmore1952) already presented a first-order approximation for the quick start of the wall movement of a bubble starting to collapse from R = Rmax after a sudden reduction of its interior pressure. Here, we present a second-order approximation for the jump start that agrees well with experimental results even under the extreme conditions of plasma-driven bubble expansion.
Models of spherical bubble dynamics including heat and mass transfer at the bubble wall (Fujikawa & Akamatsu Reference Fujikawa and Akamatsu1980; Akhatov et al. Reference Akhatov, Lindau, Topolnikov, Mettin, Vakhitova and Lauterborn2001; Lauer et al. Reference Lauer, Hu, Hickel and Adams2012; Zein, Hantke & Warnecke Reference Zein, Hantke and Warnecke2013; Peng et al. Reference Peng, Qin, Jiang and Kang2020; Zhong et al. Reference Zhong, Eshragi, Vlachos, Dabiri and Ardekani2020; Aganin & Mustafin Reference Aganin and Mustafin2021) are more complex and computationally expensive than the Gilmore model, which through its relative simplicity enables us to effectively study parameter dependencies of the bubble dynamics. Moreover, modelling of vapour condensation is hampered by the fact that the sticking coefficient for vapour molecules at the bubble wall depends on pressure and temperature, and experimental and theoretical data exhibit large variations (Marek & Straub Reference Marek and Straub2001), making predictions uncertain. Akhatov et al. (Reference Akhatov, Lindau, Topolnikov, Mettin, Vakhitova and Lauterborn2001) escaped the dilemma by using the sticking coefficient as a fit parameter for achieving a good match between the predicted rebound behaviour and experimental observations. We follow a similar strategy by using the equilibrium bubble radius after optical breakdown, Rnbd, and during collapse, Rnc 1 and Rnc 2, as fit parameters, which are chosen such that model predictions match the observed oscillation times for the first three oscillations, Tosc 1 to Tosc 3, following Vogel et al. (Reference Vogel, Busch and Parlitz1996a), Lauterborn & Kurz (Reference Lauterborn and Kurz2010) and Koch et al. (Reference Koch, Lechner, Reuter, Kohler, Mettin and Lauterborn2016). This strategy yields information on the evolution of the vapour + gas content directly after plasma formation and during bubble collapse based on a single fitting parameter for each oscillation cycle. The vapour contained in the expanded bubble at Rmax 1 is obtained by assuming equilibrium vapour pressure at room temperature, and the size the residual bubble containing mainly non-condensable gas, Rres, is derived from the frequency of the late linear oscillations. Comparison of Rnc 1 and Rres is used to quantify the vapour and non-condensable gas fractions of the collapsed bubble. Our approach cannot continuously track the evolution of gas and vapour content during bubble oscillations, as explicit models do. However, the information obtained for specific points in time enables us to assess the breakdown pressure and the collapse pressure, to simulate shock wave emission after breakdown and collapse, and to establish a complete energy balance.
The non-condensable gas in laser-induced cavitation mainly originates from water dissociation in the laser plasma, which involves free-electron-mediated processes and thermal dissociation and becomes ever more effective with increasing plasma temperature. In our experiments with tightly focused femtosecond (fs) pulses, the plasma temperature is relatively low, which leads to little permanent gas production and a very high collapse pressure.
We demonstrate the potential of our hybrid approach using experimental data on R 0 and Tosci as input for simulations with the extended Gilmore model through a detailed analysis of the scattering signal from a bubble covering more than 100 bubble oscillation cycles. The energy balance reveals that energy partitioning after optical breakdown and collapse differs strongly. During breakdown, all energy is transiently stored in the laser plasma, and a large fraction of absorbed laser energy is converted into bubble energy. By contrast, during first bubble collapse most energy goes into compression of the liquid surrounding the collapsed bubble and is then radiated away acoustically, while very little remains in the rebounding bubble. The energy of breakdown and rebound shock waves is rapidly dissipated as heat behind the shock front. This ‘convective’ heat transport involves a larger energy fraction than conductive transport by heat diffusion from the plasma or the collapsed bubble. Heating of a liquid shell around the bubble induces a local reduction of surface tension and viscosity that explains the large number of oscillations observed for highly spherical bubbles.
2. Experimental methods
The creation of highly spherical cavitation bubbles requires aberration-free, tight focusing of short laser pulses at large NA at sufficiently large distance from any solid or free surfaces to avoid jet formation (Vogel et al. Reference Vogel, Nahen, Theisen, Birngruber, Thomas and Rockwell1999a; Vogel et al. Reference Vogel, Linz, Freidank and Paltauf2008; Obreschkow et al. Reference Obreschkow, Tinguely, Dorsaz, Kobel, de Bosset and Farhat2011, Reference Obreschkow, Tinguely, Dorsaz, Kobel, de Bosset and Farhat2013; Sinibaldi et al. Reference Sinibaldi, Occhicone, Alves Pereira, Caprini, Marino, Michelotti and Casciola2019). We use long-distance water immersion objectives built into the wall of a water cell to achieve this goal. During collapse, any shape irregularities are amplified and the bubble gets deformed by a Rayleigh–Taylor instability of the bubble wall and may even disintegrate into fragments which often coalesce again during the rebound phase (Strube Reference Strube1971; Prosperetti & Hao Reference Prosperetti and Hao1999; Yuan et al. Reference Yuan, Ho, Chu and Leung2001). Observation of a large number of afterbounces without bubble disintegration is a simple practical criterion indicating that the spherical bubble shape has survived the stability crisis during the first collapse.
Besides by an elongated plasma shape of adjacent boundaries, spherical bubble dynamics may also be distorted by buoyancy if the hydrostatic pressure difference between lower and upper bubble walls and the oscillation time are large enough to induce a significant upward bubble motion during one oscillation cycle. Upon collapse, the movement is then accelerated because of the conservation of Kelvin impulse (Benjamin & Ellis Reference Benjamin and Ellis1966). This induces a fast liquid jet that penetrates the bubble and becomes visible when it rebounds after the collapse. The influence of buoyancy can be assessed by the parameter
which expresses the ratio of bubble collapse time to the time it takes an inviscid bubble of radius Rmax to rise one radius from rest driven by buoyancy forces (Blake, Taib & Doherty Reference Blake, Taib and Doherty1986; Best & Kucera Reference Best and Kucera1992). Here, ${p_\infty }$ denotes the ambient pressure, pv the vapour pressure at ambient conditions, ρ 0 is the mass density of the liquid and g the gravitational constant. The influence of buoyancy decreases with decreasing bubble size. We assume that laser-induced bubbles remain spherical up to Rmax = 100 μm (δ = 0.003), in accordance with previous observations of stable sonoluminescence for bubbles of ≈90 μm maximum radius (Gompf & Pecha Reference Gompf and Pecha2000; Weninger et al. Reference Weninger, Evans and Putterman2000). Investigation of small bubbles not only minimizes the influence of buoyancy but also provides a large dimensionless stand-off distance γ = d/Rmax from the front lens of the microscope objective. For Rmax = 100 μm, γ ≥ 22 for all focusing objectives used in this paper.
2.1. Set-up for the generation of spherical bubbles, plasma photography and oscillation tracking
The experimental arrangements for the investigation of the behaviour of spherical laser-induced cavitation bubbles in water are depicted in figure 1. We used different laser systems and detection schemes, depending on the investigation task. Small, highly spherical bubbles were generated and monitored by fs laser pulses using the set-up of figure 1(a). Their initial size was identified with the size of the laser-induced plasma and determined by photographing the plasma luminescence, whereas their oscillations were tracked by recording the scattering signal of a cw probe beam. Previous work showed excellent agreement between photographically determined R(t) curves and the predictions of spherical bubble models, including the ratio Rmax/Tosc (Kroeninger et al. Reference Kroeninger, Koehler, Kurz and Lauterborn2010; Obreschkow et al. Reference Obreschkow, Tinguely, Dorsaz, Kobel, de Bosset and Farhat2013). Therefore, it is sufficient to measure Tosc to characterize spherical bubble oscillations.
Energetic ns laser pulses were employed for time-resolved pump-probe photography of bubble wall formation and initial shock wave emission using the set-up of figure 1(c).
In figure 1(a), the laser source is a Ti:sapphire fs laser (Spectra Physics Spitfire) pumping a travelling-wave optical parametric amplifier of superfluorescence (TOPAS; Light Conversion, TOPAS 4/800) as described by Linz et al. (Reference Linz, Freidank, Liang and Vogel2016). At a wavelength of λ = 775 nm and at 1 kHz repetition rate, this laser system delivers pulses of 265 fs duration and up to 20 μJ pulse energy.
The core of the set-up for plasma photography and probe beam measurements of the subsequent bubble oscillations is a water-filled cuvette with three confocally adjusted water immersion microscope objectives (Leica HCX APO L U-V–I) built into the cuvette wall. The pump laser beam is focused into deionized and filtered (0.2 μm) water by either a $\times$40, NA = 0.8 objective with 3.3 mm working distance, or a $\times$63, NA = 0.9 objective with a working distance of 2.2 mm. The rear entrance pupil of the objectives was overfilled to create a uniform irradiance distribution corresponding to an Airy pattern in the focal plane. A cw probe laser beam (CrystaLaser, 658 nm, 40 mW) is aligned collinear and confocal with the fs-pump beam. The transmitted probe laser light is collected by a $\times$10, NA = 0.3 objective built into the opposite cell wall and imaged onto a photoreceiver (Femto HCA-S-200M-SI) connected to a digital oscilloscope (Tektronix DPO 70604). The photoreceiver is protected from the fs laser irradiation by a blocking filter.
Plasma luminescence was photographed through a third microscope objective ($\times$20, NA = 0.5, 3.5 mm working distance) that was oriented perpendicular to the optical axis of the pump and probe beams, and recorded by a digital SLR camera (Canon EOS 5D). The intermediate image formed by the $\times$20 objective and tube lens was further magnified 8 times using a Nikkor objective (63 mm/1 : 2,8). This way, we achieved a total magnification factor of 162 and a diffraction-limited spatial resolution of 0.5 μm. A confocal arrangement of all three water immersion objectives could only be achieved when the $\times$40 objective was used to focus the fs-pulses. For tighter focusing with the $\times$63 objective, the $\times$20 imaging objective had to be removed and we could only perform probe beam scattering measurements.
The absorbed fraction $E_{abs}$ of the laser energy $E_L$ was obtained from measurements of the plasma transmittance Ttra using the relation ${E_{abs}} = {E_L}(1 - {T_{tra}})$. For transmission measurements, the photoreceiver was exchanged by a calibrated energy meter, and the $\times$10 objective was replaced by a $\times$63 water immersion objective (NA = 0.9) that collected all transmitted light. Calibration accounted for light losses by reflections at optical surfaces and by absorption in the microscope objective and in water.
2.2. Single-shot recording of bubble oscillations via probe beam scattering
Figure 1(b) depicts the path of the probe laser beam through a bubble produced by a confocal pump laser pulse. The confocal adjustment guarantees that even very tiny bubbles at the optical breakdown threshold can be detected with high sensitivity (Vogel et al. Reference Vogel, Linz, Freidank and Paltauf2008; Linz et al. Reference Linz, Freidank, Liang and Vogel2016). For larger spherical bubbles, the probe beam passes perpendicularly through the bubble wall, and most of the probe laser light is transmitted through the focal region. Only a small portion is scattered or reflected at the bubble walls and interferes with the directly transmitted beam. For bubbles larger than the Rayleigh range, up to 96 % of the incident light is transmitted and 0.04 % interferes with the transmitted beam after being reflected at the rear and then the anterior bubble wall. The bias is removed by AC coupling of the photoreceiver, which had a signal bandwidth reaching from 25 kHz to 200 MHz. The coherent mixing of multiply reflected light with the transmitted beam that is shown in figure 1(b) causes a small interference modulation of the probe beam signal at the detector. Additionally, changes in the angular distribution of Mie forward scattering lead to fluctuations of the light amount transmitted through the collection aperture. These fluctuations are most pronounced for bubbles larger than the beam waist diameter but smaller than the Rayleigh range. During the oscillation of large bubbles, this leads to transient strong signal modulations shortly after the start and towards the end of each oscillation, and these modulations enable us to determine the oscillation period, Tosc (Vogel et al. Reference Vogel, Linz, Freidank and Paltauf2008). The collection NA should be chosen such that these modulations are maximized while the direct light transmission is attenuated as much as possible. In our experiments, a value of NA ≈ 0.1 provided the best results.
Due to buoyancy, small bubbles move a little upward during their oscillations although they retain a spherical shape. Nevertheless, late bubble oscillations are still detectable if the buoyancy is small enough such that the bubble stays within the probe beam focus. According to (2.1), a bubble with Rmax = 100 μm moves approximately 0.63 μm during the first oscillation but much less during later oscillations, when the radius is significantly smaller. The diffraction-limited focus diameter d = λ/NA in our set-up is 0.97 μm. Thus, the bubble moves less than the focus diameter, and late oscillations can be detected. Because the light passage through the bubble becomes slightly asymmetric after the large initial oscillation during which buoyancy is strongest, the forward Mie scattering lobe will later be obliquely oriented. Therefore, both width and orientation of the central lobe change during bubble oscillations, which enhances the intensity fluctuations of the light transmitted through the collecting aperture of NA ≈ 0.1. For 30 μm < Rmax < 50 μm, more than 100 oscillations could be traced in ‘lucky shots’, and radius variations well below 1 nm during late oscillations were detected.
2.3. Time-resolved photography of bubble formation and shock wave emission
Energetic 10 mJ and 20 mJ nanosecond laser pulses focused at NA = 0.25 were used to create large high-density plasmas, as shown in figure 1(c). Such plasmas enable us to visualize the bubble wall formation in the early phase of plasma expansion as well as shock-wave-induced phase transitions, which may enlarge the vaporized liquid volume and shift the bubble wall location. Bubbles were generated by a Nd:YAG laser (Continuum YG 671-10), which delivers pulses of 1064 nm wavelength with 6 ns duration at pulse energies of up to 250 mJ. The pulse energy was measured using a pyroelectric energy meter (Laser Precision Rj 7100).
A collimated and optically delayed frequency-doubled portion of the pump laser beam was used for illumination at short delay times up to 120 ns. For each shot, the bubble oscillation time was monitored by recording hydrophone signals of breakdown and collapse shock waves using a polyvinylidene fluoride (PVDF) hydrophone (Ceram) with 12 ns rise time. With the collimated illumination beam, we could not use the $\times$20 microscope objective for imaging as done in figure 1(a) because its back focal plane lies inside the objective and it can easily be damaged when the collimated laser beam is focused on an interior lens. Instead, we used an external $\times$7 macro objective (Leitz Photar) for photography that provided a spatial resolution of 2 μm.
3. Theoretical analysis and numerical methods
After introducing the Gilmore model of cavitation bubble dynamics, we present an extension of the model considering the rapid increase of bubble wall velocity during laser-induced energy deposition in a compressible liquid. Acoustic and shock wave emissions are then described based on the extended equation of motion. We consider water vapour generation by vaporization of the liquid in the plasma and its progressive condensation during bubble oscillations by fitting the equilibrium bubble radii such that the model predictions agree with measured values of oscillation times. Finally, we present a complete energy balance for laser-induced cavitation.
3.1. Equations governing the cavitation bubble dynamics
We used the Gilmore model of cavitation bubble dynamics (Gilmore Reference Gilmore1952; Lauterborn & Kurz Reference Lauterborn and Kurz2010) to calculate the temporal development of the bubble radius and the pressure inside the bubble, as well as the pressure distribution in the surrounding liquid. The model considers the compressibility of the liquid surrounding the bubble, viscosity and surface tension. Sound radiation into the liquid from the oscillating bubble is incorporated based on the Kirkwood–Bethe hypothesis (Cole Reference Cole1948). The Gilmore model assumes a constant gas content of the bubble, neglecting evaporation, condensation, gas diffusion through the bubble wall and heat conduction. For strong oscillations, i.e. strong compression of the contents inside the bubble, the model is augmented by a van der Waals hard core law to account for a non-compressible volume of the inert gas inside the collapsing bubble (Löfstedt, Barber & Putterman Reference Löfstedt, Barber and Putterman1993; Lauterborn & Kurz Reference Lauterborn and Kurz2010).
The bubble dynamics is described by the equation
Here, R is the bubble radius, $U = \textrm{d}R/\textrm{d}t$ is the bubble wall velocity, an overdot means differentiation with respect to time, C is the speed of sound in the liquid at the bubble wall and H is the enthalpy difference between the liquid at pressure p(R) at the bubble wall and at hydrostatic pressure
whereby ρ and p are the density and pressure within the liquid, and r is the distance from the bubble centre. The driving force for the bubble motion is expressed through the difference between the pressure within the liquid at the bubble wall and at a large distance from the wall (static pressure). Assuming an ideal gas inside the bubble, the pressure P at the bubble wall is given by
where σ denotes the surface tension, μ the dynamic shear viscosity and κ the ratio of the specific heat at constant pressure and volume. The symbol Rn denotes the equilibrium radius of the bubble at which the bubble pressure balances the hydrostatic pressure. The term $R_{v\,\textrm{d}W}^3 = {(b{R_n})^3}$ describes the size of the van der Waals hard core, with van der Waals radius $R_{v\,\textrm{d}W}$ and van der Waals coefficient b. The pressure is assumed to be uniform throughout the volume of the bubble. The pressure far away from the bubble is $p{|_{r \to \infty }} = {p_\infty }$. The equation of state (EOS) of water is approximated by the Tait equation, with B = 314 MPa, and n = 7 (Ridah Reference Ridah1988)
which leads to the following relationships for the sound velocity C and enthalpy H at the bubble wall:
with ${c_\infty }$ and ${\rho _\infty }$ denoting the sound velocity and mass density in the liquid at normal conditions. The term dH/dR in (3.1) can be derived from (3.3) and (3.6) by calculating (dH/dP) × (dP/dR). It reads
3.2. Description of laser-induced bubble initiation
Following Vogel et al. (Reference Vogel, Busch and Parlitz1996a), we neglect details of the breakdown process and refer only to the plasma size at the end of the laser pulse, and to the maximum radius reached by the cavitation bubble as a consequence of plasma expansion. Calculations start with a bubble nucleus with radius R 0, whereby the volume of this nucleus is identified with the photographically determined plasma size in the liquid. At t = 0, the nucleus contains liquid water which is then heated by the laser pulse, expands and forms a bubble when temperature and pressure have dropped below the critical point. For the sake of convenience, one usually denotes the outer border of this nucleus also already as the ‘bubble wall’. The energy input during the laser pulse is simulated by raising the value of the equilibrium radius Rn from its small initial value Rn = R 0 at the beginning of the pulse to a much larger final value Rnbd. The underlying assumption is that the absorbed laser energy is proportional to the amount of liquid vaporized by the laser pulse, which in turn is proportional to the equilibrium volume of the laser-induced bubble given by $4/3{\rm \pi} R_{nbd}^3$. This assumption holds when the energy deposited into the plasma is significantly larger than the sum of the energy needed for heating the vaporized liquid volume to the boiling temperature plus the latent heat of vaporization.
The equilibrium radius Rn is a measure of the total gas content of the cavitation bubble and does not distinguish between vapour and non-condensable gas. An increase of the Rn value beyond R 0 implies that the pressure inside the bubble rises and that the bubble starts to expand. Iteratively, we determine the Rnbd value for which the calculation yields the same oscillation time Tosc 1 or maximum radius Rmax 1 as determined experimentally.
While energy deposition by ultrashort laser pulses can be regarded as quasi-instantaneous, the finite duration of the laser pulse must be considered for ns breakdown. The temporal evolution of the laser power PL during the pulse is modelled by a sin2 function with duration τL (full-width at half-maximum) and total duration 2τL
Assuming that the cumulative volume increase of the equilibrium bubble at each time t during the laser pulse is proportional to the laser pulse energy EL absorbed up to this time, Vogel et al. (Reference Vogel, Busch and Parlitz1996a) derived an equation for the temporal development of the equilibrium radius Rn during the laser pulse
For laser-induced bubble formation in an incompressible liquid, the pressure discontinuity at the plasma border would be felt throughout the entire volume of the liquid once it is allowed to take effect in the simulation. The cavity wall then starts to be accelerated outward from rest, i.e. $U = 0$ at $t = 0$. By contrast, for a compressible liquid, the shock front represents an ‘event horizon’ up to which the breakdown effects are ‘felt’ by the liquid. As the plasma border is sharp, a shock front will form immediately, and the initial bubble wall velocity U 0 equals the initial particle velocity up behind the shock front, whereby the shock pressure is identical with the initial plasma and bubble pressure, ps = P.
Gilmore considered the case where the internal bubble pressure, Pi, is suddenly changed to a new constant value, which produces a finite velocity jump in an infinitesimal time. Considering only large terms in the equation of motion for the bubble wall and using the Tait equation that links H and C to P, he derived the first-order approximation
The approximate expression is accurate when $|H|\ll {C^2}$, which for water corresponds to $|P_{i} - {p_\infty }|\ll 2000\ \textrm{MPa}$ (Gilmore Reference Gilmore1952). This is not sufficient for modelling laser-induced breakdown, where much larger plasma pressures may be involved. Therefore, we will present a derivation of U 0 based on the Hugoniot curve data from Rice & Walsh (Reference Rice and Walsh1957). It yields (3.10) as first-order approximation and enables us to formulate a second-order approximation, which is accurate up to much higher pressure values.
Rice and Walsh fitted their Hugoniot curve data by the analytical expression
where us is the shock wave velocity and the constants are c 1 = 5190 m s−1, c 2 = 25 306 m s−1 and c 0 is the sound velocity, c∞ = 1483 m s−1. By rearranging (3.11), us can be expressed as a function of up
Using (3.11) and the conservation of momentum at a shock front, ${p_s} - {p_\infty } = {u_s}{u_p}{\rho _\infty }$ (Duvall & Fowles Reference Duvall and Fowles1963), one can link ps to us
where ${\rho _\infty } = 998\ \textrm{kg}\;{\textrm{m}^{ - 3}}$ is the mass density of water and ${p_\infty } = {10^5}\ \textrm{Pa}$ is the hydrostatic pressure. Inserting (3.12) into (3.13), one finally obtains a relation between ps and up
If ${u_p} \ll {c_1}$, the second term in the bracket can be dropped, which leads to
Immediately after breakdown, up = U 0 and ps = P. After resolving (3.15) for U 0, we get
which equals Gilmore's first-order approximation in (3.10).
Before deriving a higher-order approximation of (3.14), let us first see how we can integrate the rapid start of the bubble wall velocity during the laser pulse into the equation of motion (3.1). For this purpose, we rewrite the equation such that it describes the evolution of $\dot{U}$ and add a term ${\dot{u}_p}$ that expresses the evolution of the particle velocity at the bubble wall driven by the energy deposition during the laser pulse
The term ${\dot{u}_p}$ is derived from (3.16) as
The time interval $0 \le t \le 2{\tau _L}$ corresponds to the duration of the laser pulse as defined by (3.8). We shall now look at the pressure evolution. For very short pulse durations, energy deposition is inertially confined and we can neglect the bubble wall movement during the pulse and use the approximation R = R 0. Since the fluid does not yet move, we can neglect also viscosity. Assuming κ = 4/3, we obtain from (3.3) for the time evolution of the bubble pressure during the laser pulse
with Rn(t) given by (3.9). The time derivative of (3.19) reads
and the time derivative of (3.9) is
By inserting (3.21) into (3.20) one gets
and by inserting (3.22) into (3.18) one finally obtains
which enables us to numerically integrate (3.17).
A simulation of the bubble wall movement based on the first-order approximation is shown as dash-dotted curve in figure 2, together with simulation results not considering the ‘jump start’ of the bubble wall and the results of the second-order approximation that will be presented below. The first-order approximation yields a start velocity ${U_0} = 886\ \textrm{m}\;{\textrm{s}^{ - 1}}$, which is much higher than the particle velocity behind a shock front having a pressure equal to the bubble pressure at the end of the laser pulse. For the starting conditions of figure 2, this pressure is Pmax = 1.31 GPa, and the corresponding particle velocity is only approximately $500\ \textrm{m}\;{\textrm{s}^{ - 1}}$ (Rice & Walsh Reference Rice and Walsh1957).
Besides overestimating U 0, the first-order approximation predicts a continuous drop of the bubble wall velocity after the jump start, which contradicts the physical picture of the sequence of events. Although the shock front immediately detaches from the plasma, the bubble wall continues for a while to be accelerated by the internal bubble pressure, which leads to a peak of the U(t) curve a short while after the jump start. Later, the bubble wall velocity decreases although the bubble pressure is still higher than the hydrostatic pressure because the kinetic energy imparted to the liquid is distributed among an ever-larger liquid mass.
In order to improve the accuracy of the model predictions, we go back to the relationship between P and up in (3.14) and formulate a second-order approximation considering the second term in the bracket through its Taylor expansion
If up is well below ${c_1} = 5190\ \textrm{m}\;{\textrm{s}^{ - 1}}$, higher-order terms of the Taylor expansion can be dropped. Keeping the first term and inserting (3.24) into (3.14), we obtain
For $P \gg {p_\infty }$, we can ignore p∞, which enables us to formulate a quadratic equation of type $(A{x^2} + Bx - P = 0)$
For A, B, P > 0, such equations have a positive real and a negative imaginary root
Inserting A and B in (3.26) into the positive root of (3.27), we obtain
The time derivative of this equation is
Equation (3.29) equals the first-order approximation result in (3.18) when the second term in the denominator is neglected. Inserting (3.22) into (3.29), we finally obtain
with P given by (3.19).
Numerical integration of (3.17) with ${\dot{u}_p}$ from (3.30) yields the solid curve in figure 2, with start velocity ${U_0} = 513\ \textrm{m}\;{\textrm{s}^{ - 1}}$ in good agreement with Hugoniot data, and a time evolution U(t) that corresponds well to the expected physical scenario described above. In the following, the second-order approximation of the jump start of the bubble wall velocity will be used in all numerical simulations, if not otherwise mentioned.
3.3. Acoustic and shock wave emission
The solution of (3.17) with ${\dot{u}_p}$ from (3.30) and Rn(t) from (3.9) was used to calculate the pressure distribution in the liquid surrounding the cavitation bubble (Gilmore Reference Gilmore1952; Knapp, Daily & Hammitt Reference Knapp, Daily and Hammitt1970). The calculation is based on the Kirkwood–Bethe hypothesis, which expresses that the quantity $y = r(h + {u^2}/2)$ propagates outward along a ‘characteristic’, traced by a point moving with velocity $c + u$. Here, c is the local velocity of sound in the liquid, u is the local liquid velocity and h is the enthalpy difference between liquid at pressures p and ambient pressure p∞ (Cole Reference Cole1948). The Kirkwood–Bethe hypothesis leads to the differential equations
The pressure p at $r = r(t)$ is given by
Numerical solution of (3.17) and (3.31) with the bubble radius R, the bubble wall velocity U and the quantity $y = R(H + U{{\kern 1pt} ^2}/2)$ at the bubble wall as initial conditions yields the velocity and pressure distribution in the liquid along one characteristic. Solution of the equation for many initial conditions, i.e. along many characteristics, allows computation of u and p for a network of points (r, t). To determine u(r) and p(r) at a certain time, one has to collect a set of points with t = constant from this network.
When the bubble pressure is high, the pressure profiles in the liquid become steeper with time until a shock front is formed. Afterward, the calculations yield ambiguous pressure values, because they do not consider the energy dissipation at the shock front. The ambiguities have no physical meaning but simply indicate the presence of a discontinuity. The position of the shock front and the peak pressure at the front can be determined using the conservation laws for mass, impulse and energy flux through the discontinuity. As illustrated in supplementary figure S1 available at https://doi.org/10.1017/jfm.2022.202, it is defined by a vertical line in the u(r) plots cutting off the same area from the ambiguous part of the curve as that added below the curve (Rudenko & Soluyan Reference Rudenko and Soluyan1977; Landau & Lifschitz Reference Landau and Lifschitz1987). The location of the front was determined in the u(r) plots and transferred to the p(r) plots. The progressive reduction of peak pressure values going along with this procedure represents dissipation effects at the shock front, which are associated with an abrupt temperature rise (Brinkley & Kirkwood Reference Brinkley and Kirkwood1947; Cole Reference Cole1948; Rice & Walsh Reference Rice and Walsh1957; Duvall & Fowles Reference Duvall and Fowles1963; Müller Reference Müller2007).
We employed a commercial Matlab software package for the numerical integration of (3.17) and (3.31). The constants used for water at a temperature of 20 °C are: density of water ${\rho _\infty } = \textrm{ }998\ \textrm{kg}\;{\textrm{m}^{ - 3}}$, surface tension $\sigma = 0.073\ \textrm{N}\ {\textrm{m}^{ - 1}}$, adiabatic exponent for water vapour κ = 4/3, coefficient of the dynamic shear viscosity $\mu = 0.001\ \textrm{N s}\ {\textrm{m}^{ - 2}}$, velocity of sound ${c_\infty } = 1483\ \textrm{m}\;{\textrm{s}^{ - 1}}$, static ambient pressure ${p_\infty } = 100\ \textrm{kPa}$, vapour pressure ${p_v} = 2.33\ \textrm{kPa}$ and van der Waals coefficient b = 1/9. A van der Waals hard core is used in the calculations of bubble collapse but it is not needed for modelling the bubble expansion. Therefore, the van der Waals radius reads ${R_{v\,\textrm{d}W}} = 1/9{R_{nc}}$, where Rnc is the equilibrium radius of the bubble relevant for the collapse phase. It is considerably smaller than Rnbd immediately after optical breakdown because most of the water vapour produced during bubble generation condenses during the oscillation (Ebeling Reference Ebeling1978). The Rnc value is chosen such that the calculation yields the same oscillation time of the rebounding bubble, Tosc 2, as determined experimentally.
3.4. Choice of the adiabatic exponent
We use the room temperature value of the adiabatic exponent for water vapour, κ = 4/3, although the adiabatic exponent of water drops with increasing temperature and the temperature during breakdown and bubble collapse reaches much higher values than room temperature. This simplification is justified by the fact that part of the water molecules will dissociate for T > 3000 K (Mattsson & Desjarlais Reference Mattsson and Desjarlais2006, Reference Mattsson and Desjarlais2007), resulting in diatomic molecules such as H2 and O2, which have a larger adiabatic exponent of 1.4 at room temperature (Fujikawa & Akamatsu Reference Fujikawa and Akamatsu1980). Since both changes will, at least partly, compensate each other, the choice κ = 4/3 appears reasonable even for a large temperature range.
Bubble oscillation is isothermal (κ = 1) most of the time and adiabatic only during early expansion and late collapse/early rebound (Prosperetti & Hao Reference Prosperetti and Hao1999; Brenner, Hilgenfeldt & Lohse Reference Brenner, Hilgenfeldt and Lohse2002). Therefore, some researchers use a continuously changing time-dependent value of κ (Brenner et al. Reference Brenner, Hilgenfeldt and Lohse2002), and others switch from the isothermal to the adiabatic value, when the bubble radius passes the equilibrium radius (Barber et al. Reference Barber, Hiller, Lofstedt, Putterman and Weninger1997; Yuan et al. Reference Yuan, Ho, Chu and Leung2001). Nevertheless, following Lauterborn & Kurz (Reference Lauterborn and Kurz2010), we use a constant value because a variation of κ will largely complicate the tracking of energy partitioning during the bubble oscillations while it has little influence on R(t) and P(t).
3.5. Indirect consideration of condensation in the transition from nonlinear to linear bubble oscillations
During laser-induced plasma formation, liquid water within the plasma volume is vaporized and partially dissociated into gaseous products (Roberts et al. Reference Roberts, Cook, Rogers, Gleeson and Griffy1996; Mattsson & Desjarlais Reference Mattsson and Desjarlais2006; Elles et al. Reference Elles, Shkrob, Crowell and Bradforth2007; Müller et al. Reference Müller, Bachmann, Kroninger, Kurz and Helluy2009). Atomic hydrogen and oxygen will largely recombine to form water but some molecular hydrogen and oxygen remain as long-lived gaseous products (Nikogosyan, Oraevsky & Rupasov Reference Nikogosyan, Oraevsky and Rupasov1983; Barmina, Simakin & Shafeev Reference Barmina, Simakin and Shafeev2016, Reference Barmina, Simakin and Shafeev2017). Unlike for single bubble sonoluminescence (SBSL), where a sequence of many acoustically driven oscillations allows for rectified diffusion of dissolved air into the bubble (Brenner et al. Reference Brenner, Hilgenfeldt and Lohse2002), diffusion of dissolved gas into a laser-induced cavitation bubble is negligibly small (Akhatov et al. Reference Akhatov, Lindau, Topolnikov, Mettin, Vakhitova and Lauterborn2001). The degree of water dissociation depends on temperature and, thus, on the initial plasma energy density (Mattsson & Desjarlais Reference Mattsson and Desjarlais2006; Sato et al. Reference Sato, Tinguely, Oizumi and Farhat2013). Therefore, the vigour of the bubble collapse, which depends on the amount of non-condensable gas contained in the bubble, is correlated to the properties of the laser plasma.
The vapour produced during optical breakdown largely condenses during the nonlinear bubble oscillations, only the non-condensable gas remains, and finally the bubble exhibits small-amplitude linear oscillations around the equilibrium radius of the residual gas bubble, Rres. The use of different Rn values for the calculation of laser-induced bubble expansion and for its dynamics during the first collapse and later collapse events parametrizes vapour condensation during the first few oscillation cycles (Ebeling Reference Ebeling1978). The Rnbd, Rnc 1 and Rnc 2 values are chosen by fitting the predicted bubble dynamics to measured values of Tosc 1, Tosc 2 and Tosc 3, respectively. After the second collapse, the Rn value is kept constant because we assume that condensation is now approximately complete and that the residual bubble is mainly filled with non-condensable gas.
A reduction of Rn between breakdown and collapse seems to contradict the assumption of adiabatic expansion and collapse implied in (3.3). In fact, expansion and collapse can be approximated as adiabatic processes only during the initial expansion and the final collapse phase. In the expanded stage, heat and mass transfers at the bubble wall resemble an isothermal scenario. It is usually assumed that at R = Rmax the vapour inside the bubble is at equilibrium with the liquid outside the bubble such that the bubble pressure corresponds to the equilibrium vapour pressure at room temperature (Fujikawa & Akamatsu Reference Fujikawa and Akamatsu1980; Prosperetti & Hao Reference Prosperetti and Hao1999; Akhatov et al. Reference Akhatov, Lindau, Topolnikov, Mettin, Vakhitova and Lauterborn2001). However, the change between adiabatic and isothermal conditions hardly affects the bubble motion during the expanded stage because after the initial expansion phase, the ongoing bubble expansion is driven by inertia. Thus, condensation and heat exchange take place without any major influence on R(t) during inertially controlled oscillations but they are crucial for the final collapse phase, the collapse pressure and the rebound amplitude. Therefore, the polytropic equation (3.3) together with a reduction of Rn at the stages of maximum bubble expansion provides a realistic description of the bubble oscillation with implicit consideration of the net amount of condensation taking place during the first oscillations. The actual decrease of Rn is gradual and not stepwise as assumed in our simulations, where Rn is reduced at Rmax 1 and Rmax 2, but the stepwise reduction does not influence the predicted dynamics.
Fujikawa & Akamatsu (Reference Fujikawa and Akamatsu1980) and Yasui (Reference Yasui1995) demonstrated that the bubble collapse is more vigorous when mass transfer by condensation and heat conduction are considered in the simulations because the reduction of the bubble's gas content by condensation reduces the buffering effect of the gas. Heat conduction reduces the temperature at collapse, which facilitates condensation and leads to higher collapse velocity and peak pressure. In our approach, the influence of both condensation and heat conduction is indirectly accounted for by fitting Rnc to match measured oscillation times.
The mass reduction of the bubble content during the collapse phase must be considered for obtaining realistic values of the collapse temperature. This is done by assuming that the collapse proceeds as adiabatic process from R = Rmax starting at room temperature with a virtual bubble pressure ${p_{R\,max,virt}}$ corresponding to the amount of gas represented by Rnc. This virtual starting pressure is lower than the real pressure given by the equilibrium vapour pressure at R = Rmax. For determining ${p_{R\,max,virt}}$, we must consider that Rnc refers to a bubble with internal pressure p = p∞ at room temperature. Since ${p_{R\,max,virt}}$ also refers to room temperature, we must relate the pressure in bubbles of different size at equal temperature containing different amounts of gas, which is described by Boyle's law. That leads to
For an adiabatic collapse, pressure and temperature are linked by
which provides
for the collapse temperature. With κ = 4/3 and by inserting (3.34) into (3.36) we get
This approach provides an upper estimate of the collapse temperature, as heat conduction is neglected.
At a later stage, when it exhibits small-amplitude linear oscillations, the bubble is filled mostly with non-condensable gas. It originates largely from water dissociation in the laser plasma; Akhatov et al. (Reference Akhatov, Lindau, Topolnikov, Mettin, Vakhitova and Lauterborn2001) showed that rectified diffusion of dissolved air into the laser-induced cavitation bubble is negligibly small. Besides the non-condensable gas, a small fraction of vapour will also be present in the residual bubble. Its amount is given by the equilibrium pressure corresponding to the temperature of the liquid at the bubble wall. The linear resonance frequency of the residual bubble reads as (Lauterborn & Kurz Reference Lauterborn and Kurz2010)
Measurement of the bubble oscillation time at late stages yields ${\nu _0}$, and by inserting this value into (3.38), ${R_{nres}}$ can be determined with high precision. Comparison of the radius ${R_{nres}}$ of the residual gas bubble and Rnc then enables us to discriminate between the gas and vapour content at the first bubble collapse.
3.6. Energy balance for laser-induced bubble formation and oscillations
Already decades ago theoretical studies have shown that, during the expansion of large bubbles driven by underwater explosions or induced by optical breakdown, the largest part of the initial energy is radiated away as a shock wave and degraded into heat by dissipative processes as the wave propagates outward (Cole Reference Cole1948; Ebeling Reference Ebeling1978). A smaller fraction of the initial energy remains as bubble energy, and upon collapse and rebound of the bubble the largest part of the remaining energy is again radiated away acoustically.
Research in the 1980s and 1990s (Vogel & Lauterborn Reference Vogel and Lauterborn1988; Vogel et al. Reference Vogel, Busch and Parlitz1996a; Vogel et al. Reference Vogel, Noack, Nahen, Theisen, Busch, Parlitz, Hammer, Noojin, Rockwell and Birngruber1999b) focused on an experimental investigation of energy partitioning for millimetre-sized bubbles by measuring the bubble's potential energy
and shock wave energy
with Rm denoting the distance of the measurement location from the emission centre. Measurement of the temporal shock wave profile needed for the determination of ESW is challenging (Vogel & Lauterborn Reference Vogel and Lauterborn1988; Tinguely et al. Reference Tinguely, Obreschkow, Kobel, Dorsaz, de Bosset and Farhat2012; Lauterborn & Vogel Reference Lauterborn and Vogel2013), especially close to the source. However, both near- and far-field data are needed to assess the total emitted shock wave energy and the rate of energy dissipation upon wave propagation (Cole Reference Cole1948; Vogel et al. Reference Vogel, Busch and Parlitz1996a; Vogel et al. Reference Vogel, Noack, Nahen, Theisen, Busch, Parlitz, Hammer, Noojin, Rockwell and Birngruber1999b). One way out was to determine near-field pressure profiles by numerical calculations using the Gilmore model and far-field profiles by hydrophone measurements (Vogel et al. Reference Vogel, Busch and Parlitz1996a). Another alternative was to determine the energy dissipation from the decay of shock wave pressure with propagation distance that was obtained by measuring us(r) (Vogel et al. Reference Vogel, Noack, Nahen, Theisen, Busch, Parlitz, Hammer, Noojin, Rockwell and Birngruber1999b). These investigations revealed that up to a distance of 10 times the plasma radius, 80 %–90 % of the initial shock wave energy is dissipated.
Tinguely et al. (Reference Tinguely, Obreschkow, Kobel, Dorsaz, de Bosset and Farhat2012) established an energy balance of bubble collapse and rebound by identifying the emitted shock wave energy with the difference of bubble energies before collapse and after rebound (i.e. at Rmax 1 and Rmax 2). Considering the change of internal energy ΔUint arising from the work done by the liquid on the gas in the bubble between the two stages, the shock wave energy is ${E_{SW}} = E_{pot}^{max1} - E_{pot}^{max2} - \Delta {U_{int}}$. Experimentally, static pressure and gas pressure were varied in the ranges ${p_{stat}} \in [1,100]\ \textrm{kPa}$ and ${p_{R\,max1}} \in [1,100]\ \textrm{Pa}$, respectively. It turned out that ΔUint was negligible (<1 %) for the range of parameters investigated, which justified the approximation ${E_{SW}} \approx E_{pot}^{max1} - E_{pot}^{max2}$. The above approach is adequate for large bubbles but too simple for ${R_{max}} \to 0$, where viscous damping and surface tension must be considered. Moreover, it neglects the energy flow by water vaporization and condensation, and provides no information on the energy partitioning between shock wave emission and bubble formation after breakdown, Finally, it would be interesting to track the energy flow through the collapse phase itself, distinguishing between the energy stored in the compressed bubble content and in the liquid surrounding the bubble. In the following, we present a complete treatment of the energy flow and partitioning for laser-induced bubbles based on the Gilmore model.
3.6.1. Overview over energy partitioning
Figure 3 presents a flow diagram for the partitioning of the absorbed laser energy, Eabs. During optical breakdown, the energy absorbed in the plasma volume ${V_P} = (4/3){\rm \pi} R_0^3$ partitions into vaporization energy
and an internal energy gain ΔUint of the heated, pressurized gas volume. Here, T 1 and T 2 denote the room temperature (20°) and boiling temperature of water (100 °C), respectively, Cp = 4187 J (K kg)−1 is the isobaric heat capacity of water at 20 °C and LV = 2256 kJ kg−1 is the latent heat of vaporization at 100 °C. An equation for ΔUint will be given further below.
The expanding bubble content does work, Wgas, on the surrounding liquid, and the internal energy decreases accordingly. The index ‘gas’ refers here to both water vapour and the non-condensable gas produced by plasma-mediated water dissociation. The total energy involved in the bubble oscillation is
For isochoric energy deposition with ultrashort laser pulses, ${E_{abs}} = {E_v} + \Delta {U_{int}}$, and the work on the liquid starts only after the end of the laser pulse, when the energy of the free electrons in the laser plasma has been thermalized. In the general case, however, conversion of ΔUint into Wgas starts already during the laser pulse. During bubble expansion, the gas does work on the liquid, whereas during collapse the inrushing liquid does work on the bubble content.
During bubble expansion, the gas compresses the surrounding liquid and overcomes the liquid viscosity, the hydrostatic pressure p∞ and the pressure psurf arising from the surface tension at the bubble wall. In doing this, it creates kinetic energy Ekin of the accelerated liquid, potential energy Epot of the expanding bubble, drives the emission of a shock wave with energy ESW and does the work Wvisc by overcoming viscous damping. Altogether, the work done on the liquid involves the components
where Wstat and Wsurf denote the work done against hydrostatic pressure and surface tension.
At R = Rmax, the kinetic energy is zero and the potential energy reaches its maximum value $E_{pot}^{max1}$. At this stage, the energy of the breakdown shock wave, $E{\kern 1pt} _{SW}^{bd}$, can be obtained by evaluating all other terms in (3.42) and (3.43) and subtracting them from Eabs.
The vaporization energy Ev is needed for the phase transition itself and cannot be converted into mechanical energy of shock wave emission and bubble oscillation. It is released during the bubble oscillations by condensation of vapour at the bubble wall and heat conduction into the surrounding liquid. Besides from latent heat, the energy dissipated by condensation originates also from internal energy stored in the vapour. The total amount lost during expansion is $E_{cond}^{exp}$. The release of latent heat can be assessed by comparing the amount of vapour from the liquid in the plasma volume with the amount contained in the expanded bubble at Rmax 1.
During collapse, the potential energy of the expanded bubble partitions into energy needed to overcome liquid viscosity, a part $E_{compr}^{coll1}$ consumed for compression of the liquid surrounding the bubble and another part increasing the internal energy during the compression of the bubble content. That increase of internal energy is, however, counteracted by losses from vapour condensation. Together with the latent heat released during condensation, it constitutes the energy fraction $E_{cond}^{coll1}$. The release of latent heat is assessed by comparing the amount of vapour in the bubble at Rmax 1 with the amount contained in a bubble with equilibrium radius Rnc 1. We can distinguish between water vapour and non-condensable gas content of the bubble at first collapse by assuming that the residual bubble contains mostly non-condensable gas and that all vapour is condensed at second collapse. The gas content of the residual bubble is obtained by evaluating the late bubble oscillations with the help of (3.38), and the vapour content at first collapse is given by the volume difference Rnc 1 and Rnc 2 = Rnres.
The energy partitioning during rebound resembles the events after optical breakdown, with one crucial difference: shock wave emission is now driven not only by the re-expanding bubble content but also by the re-expansion of the compressed liquid surrounding the bubble. Therefore, a much larger energy fraction is radiated away acoustically during rebound than after optical breakdown. The parts of the rebound shock wave energy, $E{\kern 1pt} _{SW}^{reb}$, which originate from the re-expansion of the bubble and liquid are denoted ESWB and ESWL, respectively.
During later bubble oscillations (‘afterbounces’), shock wave emission turns into linear acoustic emission. We treat the energy of the acoustic waves emitted after second collapse and during later oscillations as one entity, denoted as Eacoust. Finally, a small bubble remains that contains non-condensable gas and vapour at equilibrium conditions. Since the potential energy of this bubble is zero, its energy is solely given by the residual internal energy $U_{int}^{res}$. As already mentioned, the residual bubble with radius Rnres contains mostly non-condensable gas and little vapour. The vapour content is determined by the equilibrium vapour pressure at the temperature of the liquid at the bubble wall. We will see in § 4.2.3 that this temperature is significantly larger than room temperature due to heat dissipation from the bubble content and at the front of the outgoing shock waves. However, for simplicity, we assume room temperature of the residual bubble's content when we establish the energy balance.
In the following, we provide a step-by-step account of energy partitioning.
3.6.2. Changes of internal energy and latent heat released into the liquid
For a spherical oscillating bubble containing an ideal gas under adiabatic conditions, the change of internal energy from state 1 (pgas 1, V 1) to state 2 (pgas 2, V 2) is
The energy absorbed during plasma formation is deposited into a small volume with equivalent spherical radius R 0 that is regarded as initial size of the laser-induced cavitation bubble. We assume a pressure
for the bubble nucleus at time t 0 = 0 before the laser pulse, and express the time-dependent gas pressure during the pulse through the time evolution of equilibrium radius Rn as
The internal energy of the bubble at time t during bubble expansion is then given by
Since the second term in (3.47) is very small, we will neglect it in the following.
For the bubble collapse and subsequent bubble oscillations, the van der Waals hard core must be considered. In the numerical integration of (3.17) and in the calculation of ${U_{int}}(t)$, it is introduced at the time corresponding to R = Rmax 1. After deleting the second term in (3.47) and considering the van der Waals hard core, it reads
In order to quantify the loss of internal energy of the bubble by condensation of vapour at the bubble wall during bubble expansion and collapse we first define the internal energy at the time of maximum bubble expansion for which equilibrium conditions at ambient temperature, i.e. isothermal conditions are assumed
The internal energy loss during bubble expansion is given by the difference between the energy of an adiabatically expanding bubble at Rmax (which is obtained by evaluating (3.47) at the time of maximum bubble expansion for Rn = Rnbd) and the energy corresponding to isothermal conditions at Rmax from (3.49)
In a similar fashion, the internal energy lost during bubble collapse is calculated by subtracting the energy of an adiabatically collapsing bubble with gas content corresponding to the equilibrium bubble radius at collapse (which is obtained by evaluating (3.48) for Rn = Rnc 1) from the energy corresponding to isothermal conditions at Rmax 1
The changes of internal energy during rebound and second collapse are
and
respectively. We assume that condensation is complete after second collapse such that a constant amount of residual internal energy remains.
The internal energy lost by condensation is dissipated as heat into the surrounding liquid. In addition, we need to look at the latent heat released into the liquid. For this purpose, we express the amount of vapour contained in the bubble at different instants in time (directly after the laser pulse, at Rmax 1, Rmin 1 and Rmax 2) through the radius of a vapour bubble at room temperature with pressure 0.1 MPa. The vapour bubble radius after breakdown is calculated considering conservation of mass during vaporization of the liquid in the plasma volume. With
The mass density of vapour at pv = 0.1 MPa and T = 20 °C is ${\rho _v} = 0.761\ \textrm{kg}\;{\textrm{m}^{ - 3}}$. The amount of vapour in the expanded bubble can be assessed by assuming that the vapour pressure at Rmax 1 and Rmax 2 equals the equilibrium vapour pressure at room temperature, pv = 2.33 kPa (Lauterborn & Kurz Reference Lauterborn and Kurz2010). The corresponding bubble radii for vapour at ambient pressure are then
The loss of latent heat by condensation during bubble expansion is given by
and EV from (3.41). Note, that this does not include the loss of internal energy by condensation, which will be presented in the next section. The energy transfer during the first collapse is
The calculation of $R_v^{coll1}$ is based on the assumption that at second collapse the condensation process is completed and the bubble contains only non-condensable gas. The latent heats released during rebound and second collapse are given by
and
respectively. The latent heat remaining at Rmax 2 is dissipated during the transition into linear bubble oscillations. For each stage, the total condensation loss Econd is the sum of internal energy loss and release of latent heat
3.6.3. Work done by the gas, and shock wave energies
The work done by the gas during the bubble oscillations is given by
For numerical integration, this relation must be rewritten as time integral. By transforming the differential variable from dR to dt’ using dR = (dR/dt′) × dt′ = U(t′)dt′, we obtain
where pgas is given by (3.46).
The work Wgas done on the liquid during bubble expansion (or by the liquid on the collapsing bubble) is composed of Wstat, Wsurf, and Wvisc, as expressed by (3.43). The work done at a given time to overcome the hydrostatic pressure is
The work done to overcome the surface tension is
and the potential bubble energy at any given time is therefore
Finally, the work required to overcome viscosity is
The energy balance for the conversion of the absorbed energy Eabs during the expansion process is obtained by numerical integration of the above equations up to the time tmax 1
The underscored terms are energy parts that are dissipated during expansion, the other terms remaining at t = tmax 1 are parts of the total bubble energy $E_B^{max1}$, which is the sum of the bubble's internal and potential energy. The parts of $E_{pot}^{max1}$related to hydrostatic pressure and surface tension are given by (3.59) and (3.60). Note that the energy change by condensation in (3.63) refers only to the part related to the change of internal energy; the loss of latent heat by condensation is already covered by (3.54). The breakdown shock wave energy, $E{\kern 1pt} _{SW}^{bd}$, cannot be calculated directly but it can be determined from (3.63) since all other terms are known.
The energy balance for the collapse process starts with the energy remaining at Rmin 1 and tracks their dissipation and conversion during collapse. It is obtained by numerical integration of (3.58)–(3.62) up to t = tmin 1
The underscored terms represent energy parts that are dissipated during collapse, the other terms remaining at t = tmin 1 describe parts of the total amount of energy contained in the compressed liquid and gas, $E_{compr}^{total}$. The energy stored in the compressed liquid, $E_{compr}^{coll1}$, cannot be calculated directly but it can be obtained from (3.64) because the other terms are known.
The energy balance for the rebound starts with the energy contained in the compressed bubble and liquid and tracks its dissipation and conversion during re-expansion. The balance at Rmax 2 is determined by numerical integration up to t = tmax 2
The total rebound shock wave energy, $E{\kern 1pt} _{SW}^{reb}$, is determinable from (3.69), since all other terms are known. We can even distinguish between the parts of the shock wave energy originating from the re-expansion of the compressed bubble with energy $U_{int}^{min1}$ and the compressed liquid with energy $E_{compr}^{coll1}$, which are denoted ${E_{SWB}}$ and ${E_{SWL}}$, respectively,
In order to obtain the energy balance for the afterbounces, the numerical integration is conducted up to a time at which the bubble oscillations have ceased and only a residual gas bubble remains. For this time period we have
The energy Eacoust of the acoustic radiation after second collapse and during later oscillations is calculated from the other known values in (3.67). Assuming that condensation is completed at second collapse, we have $R_{nres}^{} = R_{nc2}^{}$. Note that we need to know Rmax 3 to determine Rnc 2. For nanobubbles, experimental values of Rmax 3 may not be available. In that case, we identify the equilibrium radius during afterbounces and the residual bubble radius with Rnc 1.
The above approach for establishing an energy balance of laser-induced bubble oscillations is valid as long as heat conduction out of the energy deposition volume during the laser pulse can be neglected and when the laser pulse duration τL is much shorter than the bubble oscillation time. The characteristic thermal diffusion time for a spherical absorber is ${\tau _D} = {d^2}/8{\kern 1pt} \kappa $, where d is the focal diameter and κ is the thermal diffusivity. For pulse durations ${\tau _L} \ge {\tau _D}$, the dynamics changes from approximately adiabatic towards isothermal conditions, and the bubble expands significantly already during laser pulse. As a consequence, less work is done by the expanding gas and both acoustic radiation and the overshoot over the equilibrium radius are largely reduced. Models of such dynamics must explicitly consider heat and mass transfer. The approach presented here is valid only as long as energy deposition is thermally confined and ${\tau _L} \ll {T_{osc}}$.
4. Results
We first present plasma photographs providing R 0 data, and time-resolved photographs of the initial bubble expansion and shock wave emission. The images visualize shock-wave-induced phase transitions outside the plasma-heated region and show the process of bubble wall formation. Then one selected probe beam scattering signal is presented that traces the dynamics of a highly spherical bubble over more than 100 oscillations. This signal is analysed numerically to obtain the evolution of bubble radius, wall velocity and internal pressure, the temperature upon first collapse and the shock wave emission at breakdown and after the first bubble collapse. Examination of the transition from nonlinear to linear bubble oscillations and of late bubble oscillations provides insights into an elevated liquid temperature near the bubble wall during the late oscillations and on the relative content of vapour and non-condensable gas during the first bubble collapse. Finally, we establish a complete energy balance for the absorbed laser energy by tracking its partitioning and dissipation throughout the entire bubble lifetime.
4.1. Experiments
4.1.1. Plasma size, shock wave emission and bubble wall formation
Figure 4(a) shows luminescent plasmas in water produced by 1040 nm fs pulses of energies up to 600 nJ that were focused at NA = 0.8. The bubble threshold was at Eth,bubble = 25 nJ, corresponding to a threshold irradiance of $8.0 \times {10^{12}}\ \textrm{W}\ \textrm{c}{\textrm{m}^{ - 2}}$. For pulse energies ≥3 × Eth,bubble, luminescence could be detected photographically by integrating over many breakdown events. We identified the plasma boundary with the outer bound of the region in which luminescence could be clearly distinguished from the uniform background. The volume of the luminescent region determined from the photographs is plotted in figure 4(b) as a function of laser pulse energy, and figure 4(c) shows the corresponding radius values of spheres with same volume that define the initial bubble radius R 0 for the numerical simulations. Figure 4(d) shows the pulse energy dependence of plasma transmittance Topt, from which the absorbed energy is determined as ${E_{abs}} = {E_L}(1 - {T_{tra}})$.
In the theoretical description of the bubble dynamics, we assume a homogeneous pressure distribution inside the laser-produced bubble throughout the entire bubble lifetime. The initial bubble wall position is identified with the outer boundary of the luminescent plasma region, and it is assumed that the location of the bubble wall is affected only by the pressure difference between inner and outer pressure but not shifted by phase transitions arising from the shock wave passage. Both assumptions are checked by evaluating the time-resolved photographs of the initial hydrodynamic processes shown in figure 5.
First, we determined the plasma energy density, ε, by relating the luminescent plasma volume determined from the photographs to the amount of energy absorbed in this volume that is obtained from transmission measurements (Nahen & Vogel 1996). The average energy density is $\varepsilon \approx 40\ \textrm{kJ}\ \textrm{c}{\textrm{m}^{ - 3}}$ for the 10 mJ pulse in figure 5(a), and $\varepsilon \approx 35\ \textrm{kJ}\ \textrm{c}{\textrm{m}^{ - 3}}$ for the 20 mJ pulse in figure 5(b). Under the assumption of isochoric energy deposition, we can derive the plasma pressure from the energy density and obtain values of 11.3 GPa and 10.1 GPa, respectively, using the IAPWS-95 formulation of the water EOS (Wagner & Pruss Reference Wagner and Pruss2002). These data are an upper estimate; the actual pressure is somewhat lower because the bubble starts to expand already during the ns laser pulse. The large pressure jump at a shock front results in rapid energy dissipation and in a temperature rise after shock wave passage, which for Δp = 10 GPa amounts to 576 °C (Rice & Walsh Reference Rice and Walsh1957). Therefore, the dissipated shock wave energy will create a phase transition in a thin zone extending beyond the rim of the expanding plasma, up to which vaporization is induced directly by the absorbed laser energy. Based on this background information, we will now step-by-step analyse the image series.
Plasma luminescence is visible on all images although it rapidly ceases after the end of the pulse because the photographs were taken with open shutter in a darkened room. The time given on the images refers to the time delay between the pump pulse producing the plasma and the illumination pulse for shadowgraph photography. Photographic exposures were adjusted to provide similar background brightness in all pictures. However, due to the divergence of the illuminating laser beam, the illuminated spot is larger for longer delay times. This lowers its irradiance and makes the self-luminescent plasma appear brighter at late times although its actual luminescence remains the same.
Plasma formation starts at the laser beam waist (on the left) and moves upstream towards the incoming laser beam while the laser power increases during the pulse (Docchio et al. Reference Docchio, Regondi, Capon and Mellerio1988; Vogel et al. Reference Vogel, Nahen, Theisen and Noack1996b). The movement of the breakdown wave results in a delayed shock wave emission from the upstream part of the plasma. The geometrical form of the shock wave reflects both the overall plasma shape and inhomogeneities of the energy density distribution within the plasma. Such inhomogeneities can result in the release of pressure transients propagating in the zone between plasma and outer shock front that are visible as dark structures on the images. Because of the pressure dependence of sound velocity, the transients propagate at high speed and finally catch up with the shock front.
In figure 5(a), the plasma exhibits an inhomogeneous energy distribution with a high-density region close to the beam waist and a larger and more pronounced hot spot located further upstream. The inhomogeneity is due to the temporal pulse shape of the 6 ns pulse exhibiting two peaks (Vogel et al. Reference Vogel, Nahen, Theisen and Noack1996b). Therefore, the breakdown wave moving upstream against the incoming laser beam produces two spatially separated high-density regions corresponding to the two peaks of the laser pulse. The hot spot in the upstream region is generated during the second half on the laser pulse, and the fast transient emerging from it propagates into the high-pressure region behind the shock wave emitted earlier during the pulse near the beam waist. The velocity of this transient is indicative of the speed of pressure equilibration within the breakdown region and, due to the pressure dependence of sound speed, it provides information about the average pressure in this region. The pressure transient traversing the breakdown region in axial direction has passed the plasma region after approximately 40 ns, and its average speed in axial direction during the first 60 ns is approximately 300 μm/60 ns, i.e. ${\approx} 5000\ \textrm{m}\;{\textrm{s}^{ - 1}}$. For the water Hugoniot centred at ambient conditions (20 °C and 0.1 MPa), this value corresponds to a pressure of ≈10 GPa (Rice & Walsh Reference Rice and Walsh1957), which is consistent with the pressure value derived from the plasma energy density. It is also consistent with pressure values obtained by analysis of the initial shock wave speed close to the plasma rim through time-resolved photography (Vogel et al. Reference Vogel, Busch and Parlitz1996a) and streak photography (Noack & Vogel Reference Noack and Vogel1998). The shock passage results in rapid pressure equilibration within the breakdown region (a passage time of 60 ns corresponds to the 1/2800 part of the bubble expansion time of 168 μs, which was obtained from the hydrophone signal). This justifies the assumption of a homogenous bubble pressure made in the Gilmore model.
A black region between the luminescent plasma region and the shock wave appears at t = 7 ns, shortly after the peak of the 6 ns laser pulse. The outer border of this region is usually identified with the wall of the expanding cavitation bubble but it becomes visible already before a phase boundary is formed in the fluid. A phase change (i.e. the formation of the bubble wall) occurs when and where pressure and temperature both drop below the critical point (374 °C, 22 MPa). For the 10 mJ, 6 ns laser pulse, this happens only after 70–80 ns, as shown by numerical simulations of p(t) in Vogel et al. (Reference Vogel, Busch and Parlitz1996a). What we see as the ‘bubble wall’ before the formation of a phase boundary is actually a mass density jump between the hot and rapidly expanding supercritical fluid from the plasma region and the surrounding colder liquid compressed by the shock wave passage. Since lowering of the mass density goes along with a reduction of the refractive index, the optical properties of the supercritical fluid region resemble those of a bubble. For the sake of simplicity, the sharp grey level transition on the shadowgraphs is taken as the ‘bubble wall’ already from the time of laser exposure although a phase boundary appears only later.
Due to the energy dissipation at the shock front, a thin layer around high-density plasmas is heated to temperatures above the boiling point or even the superheat limit. The heated region shows up already shortly after the shock wave passage. It is clearly visible in figure 5(b), where the plasma was produced with a 20 mJ laser pulse. Here, the outer border of the heated region appears rugged, especially in the picture taken at t = 12 ns. Where the temperature rises above the superheat limit (kinetic spinodal), which lies at ≈300 °C (Debenedetti Reference Debenedetti1996; Vogel & Venugopalan Reference Vogel and Venugopalan2003; Vogel et al. Reference Vogel, Noack, Huettman and Paltauf2005), the instable liquid starts to expand immediately (Zhigilei et al. Reference Zhigilei, Leveugle, Garrison, Yingling and Zeifman2003). Therefore, the refractive index drops, which become visible as a dark region on the photographs, similar to the supercritical fluid in the plasma region. The rugged appearance of the borderline of the shock-induced phase transition is likely due to inhomogeneities in the plasma producing pressure fluctuations in the near field and local variations of the temperature rise. Thus, the shock wave can enlarge the volume of vaporized liquid driving the bubble expansion. This ‘convective’ heat transport by shock wave propagation and energy dissipation at the shock front is much faster than heat diffusion. It produces a thermal gradient that is less steep than the initial gradient arising from the spatial distribution of free electron density in the laser plasma. As a consequence, the ‘bubble wall’ appears fuzzy in figure 5 during the first 50–60 ns. It smooths out later during the ongoing plasma expansion, when a phase boundary is formed, the bubble content cools down adiabatically and surface tension can smooth out local irregularities.
The additional vapour mass produced by energy dissipation at the shock front is hard to quantify and, therefore, not included in the energy balance presented in this paper. It will have little influence on bubble expansion, which is driven by the plasma pressure that by far exceeds the pressure evolving through shock-wave-induced heating. However, the temperature increase in the liquid around the bubble reduces viscosity, which lowers the damping during its late oscillations (see § 4.2.3 further below).
4.1.2. Single-shot probe beam signal
Figure 6 shows a probe beam signal representing the oscillations of an almost perfectly spherical bubble with 35.8 μm maximum radius, a dimensionless stand-off distance γ = 70 from the microscope objective's front lens and very little influence of buoyancy (δ = 0.0017). The signal covers 102 oscillations and portrays the transition from initial nonlinear cavitation bubble oscillations to linear oscillations around the equilibrium radius of the residual gas bubble.
The bubble was produced by a 755 nm, 155 nJ fs pulse focused at NA = 0.9. The equivalent spherical plasma radius for this pulse energy read from figure 4(c) is R 0 = 1.33 μm, and the absorbed energy is 81 nJ. With the plasma volume ${V_p} = 9.85\;\mathrm{\mu }{\textrm{m}^3}$ from figure 4(b), this yields an average volumetric energy density of 8.73 kJ cm−3. The internal energy density ${U_{int}}/{V_p} = ({E_{abs}}-{E_v})/{V_p}$ with Ev = 25.49 nJ from (3.41) is 6.14 kJ cm−3, and the difference of both values corresponds to the vaporization enthalpy. Since energy deposition is isochoric, the average plasma temperature can be determined from Uint/Vp and the water EOS. Using the IAPWS-95 formulation (Wagner & Pruss Reference Wagner and Pruss2002), we obtain Tavg = 1550 K. The peak temperature will be somewhat larger, as suggested by the inhomogeneous brightness distribution in the photographs of figure 4(a).
The slightly elongated plasma shape causes a stability crisis in the first collapse that is reflected by the asymmetry of the probe beam signal during second and third oscillations seen in figure 6(b). However, the signal symmetry is regained in the fourth oscillation, which indicates that surface tension has restored the spherical shape.
The signal undulations visible during the first cavitation bubble oscillation are interference fringes reflecting the radius–time evolution. However, the region with detectable fringe separation is too small to gain significant information about the R(t) curve. Later, each signal undulation represents one period of the small-amplitude oscillations of the residual bubble around its equilibrium radius. The undulations arise from the interference between bubble wall reflections with the directly transmitted beam combined with changes in the angular distribution and orientation of the central Mie scattering lobe, as described in § 2.2. The slow undulation of the average signal level within the first 40 μs is a consequence of the capacitive AC coupling having a lower cut off frequency at 25 kHz and has no physical meaning. It does not affect the determination of bubble oscillation times. Oscillation times rapidly drop from Tosc 1 = 6.5 μs to values below 1 μs and converge against a value of 813 ± 6.7 ns during the late gas bubble oscillations.
4.2. Numerical calculations
4.2.1. Evolution of bubble radius, wall velocity and pressure; collapse temperature
Table 1 summarizes characteristic breakdown and bubble parameters corresponding to the probe beam signal of figure 6. Figure 7 shows simulation results for the time evolution of cavitation bubble radius, internal pressure and bubble wall velocity obtained with these parameters. Enlarged views of R(t), P(t) and $U(t)$ for time intervals of 20 ns after optical breakdown and 1 ns around the first collapse are presented in figure 8.
Because of the relatively small plasma temperature of 1550 K, the breakdown pressure is merely 1.25 GPa, much lower than in previous experiments with bubbles generated by 10 mJ IR ns laser pulses (Vogel et al. Reference Vogel, Busch and Parlitz1996a). Water dissociation relies on free-electron-mediated pathways as thermal dissociation sets in only at 3000 K (Mattsson & Desjarlais Reference Mattsson and Desjarlais2006). Compared with the amount of vaporized liquid, only a small amount of non-condensable gas is produced by dissociation. As a consequence, the bubble collapse is only weakly damped by its permanent gas content, and the collapse pressure reaches 13.5 GPa, which is approximately 11 times larger than the plasma pressure upon breakdown. The minimum bubble radius at collapse is Rmin = 419 nm, less than one third of the plasma size that defines the initial bubble radius. The duration of the collapse pressure peak is much shorter with a full width at half maximum (FWHM) of 76.5 ps than the pressure peak after breakdown (FWHM = 920 ns). It is interesting to note that models of bubble collapse in compressible liquids that are based on full solutions of the Navier–Stokes equations rather than on the Kirkwood--Bethe hypothesis yield somewhat higher collapse pressures than the Gilmore model (Fuster, Dopazo & Hauke Reference Fuster, Dopazo and Hauke2011; Koch et al. Reference Koch, Lechner, Reuter, Kohler, Mettin and Lauterborn2016). Thus, the actual collapse pressure may be even larger than 13.5 GPa.
During breakdown, the bubble wall velocity performs a jump start and reaches a peak value of 540 m s−1 (figure 8e). At this time, the bubble content is still a supercritical fluid. A phase boundary between vapour and liquid water forms after 6.85 ns, when pressure and temperature inside the bubble drop below the critical point (figure 8c). The bubble expansion velocity reaches its peak already after 500 ps, when the internal pressure of the plasma/bubble region still imparts kinetic energy to the surrounding liquid. Nevertheless, the velocity of the boundary between both regions already decreases because an ever-larger mass is involved in the radial outward flow. The flow stops, when at R = Rmax all kinetic energy has been converted into potential energy of the expanded bubble.
During collapse, the bubble wall is accelerated and its velocity becomes supersonic with respect to the sound velocity in the gaseous bubble content at ambient conditions, c 0. This happens at a bubble pressure of 0.45 MPa, when the mass density of the bubble content is still relatively low. The growing bubble pressure in the final collapse phase rapidly reverses the direction of the bubble wall velocity. It changes within ≈40 ps from a peak collapse value of −1788 m s−1 to a peak rebound value of 369 m s−1 (figure 8f) which goes along with an acceleration of 5.4 × 1013 m s−2. During the final collapse stage, the bubble content becomes supercritical, the phase boundary at the bubble wall disappears and reappears again during rebound.
When the bubble content becomes a supercritical fluid and is compressed to a state resembling the van der Waals hard core, the sound velocity inside the bubble rapidly increases to value larger than the bubble wall velocity. The sound velocity can be estimated by applying thermodynamic data for compressed liquid water such as presented by Rice & Walsh (Reference Rice and Walsh1957) to the compressed bubble content. Unfortunately, their data for the water Hugoniot centred at 20 °C and 0.1 MPa provide only a rough estimate because the temperature in the collapsed bubble is much higher than for the 20 °C Hugoniot at a pressure of 13.5 GPa. However, the image series in figure 5(a) provides additional information because it indicates that the pressure transient emitted from the hot spot propagates even faster through the hot plasma region than through the cooler liquid around. Thus, c will be of the order of 5000 m s−1 or faster also in the collapsed bubble. The high sound velocity promotes a rapid pressure equilibration, which again justifies the assumption of a homogeneous bubble pressure made in the Gilmore model.
The collapse temperature calculated using (3.37) for an adiabatic collapse with an amount of vapour corresponding to Rnc 1 is Tcoll = 31 400 K. This is more than ten times larger than what is found under the assumption that the collapse starts with the equilibrium vapour pressure at room temperature and proceeds with constant vapour content, neglecting condensation. Under those conditions, the buffering by the larger vapour content results in a collapse pressure pcoll = 0.061 GPa, and a temperature of 2995 K.
4.2.2. Shock wave emission at breakdown and upon bubble collapse
Figures 9 and 10 present the evolution of the velocity distributions u(r) and pressure distributions p(r) after the optical breakdown and during the final collapse and early rebound phase, respectively. After breakdown, a shock front forms within a few picoseconds, owing to the ultrashort laser pulse duration. The shock wave detaches within ≈10 ns, when a velocity and pressure minimum have evolved between the shock front and the bubble wall region. The situation is more complex in the case of bubble collapse and rebound, as seen in figure 10. The velocity of the inrushing flow reaches a maximum ≈75 ps before collapse. The high bubble pressure first stops the inward flow at the bubble wall and then drives an outward flow that reaches its peak velocity 90 ps after the collapse at a location slightly ahead of the bubble wall. The outward flow collides with the still incoming flow from outer liquid regions at the location of the shock front. Thus, the change of motion of the bubble wall is communicated to the liquid by the passage of the shock wave, and the flow reversal occurs at the shock front.
For the transients emitted after breakdown and upon the bubble's rebound, the pressure decay is faster than for acoustic waves for which the attenuation would be proportional to r −1. The steepest slope is −1.27 for the breakdown shock wave, and −1.75 during rebound. This is because the collapse pressure is one order of magnitude larger than the plasma pressure in fs breakdown. In both cases, the shock front persists also in the far field, where the pressure has dropped to a few MPa. Here, the peak amplitude decays proportional to approximately r −1.15, which is typical for the regime of weak shock wave propagation (Arons Reference Arons1954; Rogers Reference Rogers1977).
The dissipation rate of the shock wave energy is largest close to the plasma or collapsed bubble because it is proportional to the pressure jump at the shock front (Vogel et al. Reference Vogel, Noack, Nahen, Theisen, Busch, Parlitz, Hammer, Noojin, Rockwell and Birngruber1999b). Upon rebound, a shock front develops within 50 ps and 750 nm propagation distance. At this time, it exhibits a pressure jump of ≈8 GPa, corresponding to a temperature jump to 436 °C (Rice & Walsh Reference Rice and Walsh1957). Most of the shock wave energy is dissipated during the first 2 μm propagation distance, where the pressure drops to 1/10 of its peak value and the pressure decay curve is steepest (figure 10b). This heats a few micrometre thick liquid layer around the expanding bubble, which vaporizes a thin liquid shell and reduces surface tension and local viscosity in the liquid near the bubble wall.
4.2.3. Transition from nonlinear to linear bubble oscillations
Table 1 summarizes bubble parameters for the simulations of figures 7–10. We see that during the initial nonlinear oscillations most of the vapour produced during plasma formation condenses. At second rebound, the vapour content has already dropped to 1/200 of the initial value. This finding justifies our assumption, that condensation is complete at second collapse and the residual bubble undergoing linear oscillations contains only non-condensable gas.
It is interesting to note that Rv after plasma formation is larger than Rnbd. The entity Rnbd does not express the exact amount of gas and vapour in the bubble but rather measures the bubble's internal energy that does work on the surrounding liquid (see (3.48)). The difference between Rv and Rnbd is small for high-density plasmas but it becomes ever larger when the pulse energy is reduced and the plasma energy density decreases. Upon first collapse, the vapour content of the bubble is already much smaller than after breakdown, and the amount of non-condensable gas becomes relevant for buffering the collapse. Nevertheless, we see that Rvc 1 is not much smaller than Rnc 1, which reflects the total gas content of the bubble. This indicates that also water vapour significantly contributes to buffering the collapse.
We will now use the Rnc data from table 1 to estimate the ratio of vapour and permanent gas at first collapse. During the first collapse, some vapour remains in the bubble because in the final collapse stage, condensation cannot keep up with the rapid reduction of the surface area (Storey & Szeri Reference Storey and Szeri2000). However, during second collapse almost all vapour will condense because the rebounded bubble at Rmax 2 contains much less vapour than the larger bubble at Rmax 1. Thus, we can assume that the bubble content at second collapse is almost exclusively non-condensable gas. Its amount is given by the equilibrium radius Rnc 2 = 2.44 μm, while Rnc 1 = 3.6 μm stands for a gas-vapour mixture. Comparison of both values yields a vapour/gas ratio of 68.9 % vapour to 31.1 % non-condensable gas for the first collapse.
The radius of the residual bubble, Rres, can be derived from the frequency of the bubble oscillations at late times using (3.38). The mean oscillation time from 50th to 100th oscillation is 813.3 ± 6.7 ns, which corresponds to an oscillation frequency of 1.23 MHz. With room temperature values of surface tension and viscosity this yields a value Rnres = 2.60 μm for the equilibrium radius of the residual gas bubble under isothermal conditions (κ = 1). Note that Rnres = 2.60 μm, is slightly larger than Rnc 2 = 2.44 μm. We assume that during the second collapse almost all vapour remaining after the first collapse condenses at the bubble wall, whereas later some re-evaporation into the residual bubble occurs.
Figure 11 presents simulation results for the R(t) curve up to the end of the experimentally observed time period assuming different bubble wall temperatures, TW, and either adiabatic or isothermal conditions. Thermal dissipation by heat conduction from the bubble and energy dissipation from the shock wave raises the temperature of the liquid around the oscillating residual bubble. The heated liquid shell around the bubble is relatively thick during the late oscillations when the bubble is small, and the conditions are, thus, isothermal. By contrast, the heated boundary layer is very thin during the first oscillations, when the bubble is large, and the bubble dynamics can here be described well under the assumption that the bubble oscillates in water at room temperature.
The simulation results for room temperature (TW = 20 °C) predict much stronger damping than observed experimentally (figure 11a). Damping is weakest for the isothermal curve but even here, the peak-to-peak oscillation amplitude drops below 1 nm after little more than 30 μs (30 oscillations). Figures 11(b) and 11(c) show R(t) curves for bubble wall temperatures of 60 °C and 110 °C, where surface tension and, particularly, viscosity are much lower than at room temperature (see supplementary figure S2). It turns out that an average bubble wall temperature TW = 110 °C provides good agreement with experimental results.
An upper bound for the possible bubble wall temperature is given by the condition that the equilibrium vapour pressure must be balanced by the sum of hydrostatic and Laplace pressure from surface tension; otherwise, the bubble would grow in size. At Rres = 2.6 μm, this is the case for T = 113 °C, slightly above the value assumed in figure 11(c).
The possible temperature rise can also be estimated from the absorbed laser energy Eabs = 62.5 nJ (see § 4.2.4 below). It can produce an average temperature rise of 48.7 K in a 2 μm thick liquid shell around a bubble with 2.5 μm radius, corresponding to a final temperature of 68.7 °C. The thermal relaxation time for a spherical source of 9 μm diameter is 71 μs, which is similar to the observed oscillation time. Since the bubble wall temperature is larger than the average temperature in the liquid shell, the estimate Tavg = 68.7 °C is consistent with a value TW = 110 °C providing good agreement with experimental results.
The above estimate provides only a rough assessment of the influence of the elevated temperature in the bubble's surroundings as they neglect heat diffusion. Nevertheless, they show that a bubble wall temperature >100 °C is needed to explain experimental R(t) data. The predicted oscillation amplitude at Tosc, 100 is then below 0.1 nm, just detectible by our interferometric probe beam technique.
4.2.4. Energy partitioning
Table 2 presents the energy flow during bubble expansion, collapse, rebound and afterbounces for the signal of figure 6. In each column, first the energy losses are listed that in figure 3 have been depicted as dashed lines, followed by a categorization of the energy remaining at the end of the respective oscillation phase (marked by solid lines in figure 3). Only at these instants (Rmax 1, Rmin 1 and Rmax 2), it is possible to establish a complete energy balance because the fluid movement has stopped, Ekin = 0, and the amount of gas and vapour can be assessed based on simple assumptions without explicit modelling the kinetics of heat and mass transfer during the bubble oscillations. Figure 12 shows the time evolution of those energy fractions that can be continuously tracked, with emphasis on the partitioning of internal bubble energy, i.e. that part of the deposited energy which can do work on the surrounding liquid.
The incident energy of the laser pulse producing the bubble with Rmax = 35.8 μm was 155 nJ, and the absorbed energy determined from the measured plasma transmission was Eabs = 81 nJ (§ 4.1). The total amount of deposited energy derived from the plasma radius R 0 taken from figure 4 and the model fit of R(t) to measured oscillation times is ${E_{abs}} = {E_v} + \Delta {U_{int}} = 25.5\ \textrm{nJ} + 37.0\ \textrm{nJ} = 62.5\ \textrm{nJ}$. The discrepancy can be explained by the fact that the measurement accounts for the light transmitted into the aperture of the objective in front of the photodetector but cannot distinguish between absorption and scattering losses. Thus, the simple calculation ${E_{abs}} = {E_L}(1 - {T_{tra}})$ provides an upper limit for the absorbed laser energy.
Because of the moderate plasma energy density, 40.8 % of the absorbed energy is needed to vaporize the liquid in the plasma volume, and only 59.2 % are available for doing work on the surrounding liquid. A fraction of 55.4 % of $\Delta {U_{int}}$ (of Eabs) appears as potential energy $E_{pot}^{max1}$ of the expanded cavitation bubble at Rmax 1, and 38.9 % of $\Delta {U_{int}}$ (of Eabs) are emitted as shock wave. The splitting ratio of absorbed laser energy into bubble and shock wave energy is 32.8 % versus 23.0 %, i.e. 1.43 : 1.
The picture is very different for the rebound phase, where only 1.94 % of the energy goes into bubble expansion and 96.3 % into shock wave emission. The difference can be understood by looking at the energy stored upon collapse in the compressed bubble content and the surrounding liquid. It turns out that the energy content of the compressed liquid, $E_{compr}^{coll}$, is 12.6 times higher than the inner energy of the bubble at Rmin. Therefore, most energy is radiated away acoustically upon rebound and only a small fraction originating from the internal energy of the compressed bubble content contributes to bubble formation. As a consequence, the fraction of the rebound shock wave energy $E{\kern 1pt} _{SW}^{reb}$ originating from the re-expansion of the compressed liquid $({E_{SWL}})$ is 18.45 times larger than the part provided by the rebounding bubble $({E_{SWB}})$.
The energy consumed for vaporization of the liquid in the plasma volume is during the bubble oscillation dissipated via condensation and heat conduction. From the internal energy, 91.5 % is carried away by shock waves emitted after optical breakdown and bubble collapse, and only 4.7 % are lost by viscous damping. A small fraction of $\Delta {U_{int}}$ (3.8 %) is dissipated by condensation of the vapour inside the bubble. This adds to the energy required for vaporization of the liquid in the plasma volume that is later released by condensation. Although a large fraction of the deposited laser energy appears transiently as mechanical energy of the cavitation bubble and the breakdown and collapse shock waves, most energy is soon transformed into heat. We have seen in § 4.2.2 that even the shock wave energy is not just radiated away but much of it is dissipated as heat in close vicinity to the optical breakdown or collapse site. Note that in laser surgery, the bubble and shock wave energy will partly do mechanical work on cellular and tissue structures instead of being thermally dissipated.
5. Discussion
The detailed characterization of spherical cavitation bubble dynamics in this paper relies on a ‘hybrid’ approach in which numerical simulations with an extended Gilmore model are fitted to experimental data on plasma size and bubble oscillation times. In § 5.1, we first compare our hybrid approach with previously published explicit models of bubble generation, and in § 5.2 we discuss the origin of non-condensable gas in the bubble by plasma-mediated water dissociation. In § 5.3, we then describe how the amount of non-condensable gas produced in the laser plasma influences the vigour of the bubble collapse and compare the extreme conditions produced by the collapse of laser-induced bubbles with the conditions created during acoustically driven SBSL. When the laser-induced bubble contains little non-condensable gas, it collapses very violently and emits shock waves exhibiting a much more rapid pressure decay than acoustic waves. In § 5.4, this result is compared with the findings in previous studies on acoustic emission by collapsing and rebounding bubbles. Our approach can distinguish between energy fractions stored in the compressed content of the collapsed bubble and the compressed liquid surrounding the bubble, and track their conversion into bubble and shock wave energy during rebound. In § 5.5, we use this feature to elucidate the differences between energy partitioning after optical breakdown and bubble collapse, and compare our findings with the results of explicit approaches based on the solution of the Navier--Stokes equations.
5.1. Bubble generation
In ‘ab initio’ models of laser-induced cavitation, the modelling starts by describing the process of nonlinear energy deposition and subsequent phase transitions leading to bubble formation and continues with the bubble dynamics driven by the phase transition (Glinsky et al. Reference Glinsky, Bailey, London, Amendt, Rubenchik and Strauss2001; Byun & Kwak Reference Byun and Kwak2004; Dagallier et al. Reference Dagallier, Boulais, Boutopoulos, Lachaine and Meunier2017). In the present paper, the process of bubble generation is not explicitly modelled. Instead, the expansion of a virtual seed bubble is linked to two experimental parameters: the breakdown volume derived from photographs of plasma luminescence (or from the size of the breakdown region in shadowgraph photos), and the duration of the bubble oscillations determined from probe beam scattering. The seed bubble radius R 0 is obtained by identifying its volume with the breakdown volume, and model predictions for Tosci are fitted to experimental values by varying the equilibrium bubble radii Rnbd, Rnc 1 and Rnc 2, which represent the driving forces for the first, second and third oscillation, respectively. This way, comprehensive information about the laser-induced bubble dynamics becomes available; including the amount of absorbed laser energy and its partitioning, phase transitions during bubble generation and heat and mass transfers during the bubble oscillations. The latter information is obtained through the fitting procedure, without explicit modelling of these processes.
While ab initio models directly describe the evolution of the pressure driving the bubble expansion, the pressure evolution is in our model encoded in the temporal evolution of the equilibrium radius Rn. This approach can describe bubble formation in a large range of laser pulse durations, whereas explicit modelling becomes difficult for ps and ns pulses, where plasma formation, phase transition and bubble expansion occur simultaneously during the laser pulse.
While even simple, incompressible models provide an excellent description of the inertial bubble oscillation, liquid compressibility must be considered to adequately describe the interplay between bubble wall movement and shock wave emission. We extended the equation of motion of the Gilmore model by a term describing the initial jump of the bubble wall velocity to a large finite value equal to the particle velocity behind the plasma-induced shock front. The jump start of the bubble wall was observed experimentally by Vogel et al. (Reference Vogel, Busch and Parlitz1996a) but at that time not yet included in the modelling of bubble generation. Figure 13 compares predicted U(t) curves with and without jump start to experimental data from that paper. For all investigated laser parameters, consideration of the jump start provides much better agreement between experimental data and simulated U(t) curves. It is interesting to note that for the bubble in figure 13(d), where a 10 mJ, 6 ns pulse produced an average plasma energy density of 40 kJ cm−3, both measured and predicted bubble wall velocities exceed the sound velocity in the liquid. In this case, the Keller–Miksis model that considers compressibility assuming a constant sound velocity in the liquid would not be able to describe the laser-induced bubble expansion and shock wave emission correctly.
Measured peak velocities are larger than predicted values, even when the jump start is considered. This is because laser pulses were only moderately focused at NA = 0.25, which resulted in elongated plasmas, especially at higher pulse energies. Under these conditions, the bubble wall velocity in the direction perpendicular to the long plasma axis is initially faster than for spherical symmetry. In a later stage of the bubble expansion, when the bubble acquires an approximately spherical shape, the bubble wall velocity is equal in all directions. In the transition phase, the measured velocity must be slower than for spherical expansion to compensate for the initial overshoot. This behaviour is indeed observed in all cases presented in figure 13: experimental U(t) data initially exceed the predicted values for spherical expansion, then drop below the simulated U(t) curve and finally both curves converge.
Figure 14 compares the evolution of bubble pressure and shock wave profiles with and without bubble wall jump start for the 10 mJ, 6 ns pulse of figure 13(d). The shock front forms within 8 ns in both cases but the jump start of the bubble wall results in a considerably lower maximum bubble pressure (4.75 GPa) compared with the value obtained without consideration of the particle velocity behind the shock front (8.8 GPa). The experimental value for the pressure at the plasma rim was 7.15 GPa. It was obtained from the initial shock front velocity of 4500 m s−1 using (3.11) that is based on Hugoniot data valid up to 25 GPa (Rice & Walsh Reference Rice and Walsh1957). The Tait equation (3.4) used in the Gilmore model leads to the relationship
between velocity and pressure at the shock front (Müller Reference Müller1987; Vogel et al. Reference Vogel, Busch and Parlitz1996a). For us = 4500 m s−1, (3.11) yields a 1.63 times larger pressure than (5.1). If we correct the Gilmore model predictions by that factor, we obtain ps = 7.73 GPa, in very good agreement with the experimental value of 7.15 GPa
5.2. Evolution of vapour and non-condensable gas content
The strength of spherical bubble collapse depends on the amount of gas and water vapour that buffers the inrushing liquid motion (Fujikawa & Akamatsu Reference Fujikawa and Akamatsu1980; Storey & Szeri Reference Storey and Szeri2000; Akhatov et al. Reference Akhatov, Lindau, Topolnikov, Mettin, Vakhitova and Lauterborn2001; Zein et al. Reference Zein, Hantke and Warnecke2013; Zhong et al. Reference Zhong, Eshragi, Vlachos, Dabiri and Ardekani2020; Trummler, Schmidt & Adams Reference Trummler, Schmidt and Adams2021). The specific nature of the non-condensable gas will differ, depending on the type of cavitation. Akhatov et al. (Reference Akhatov, Lindau, Topolnikov, Mettin, Vakhitova and Lauterborn2001) showed that diffusion of dissolved gas into a laser-induced cavitation bubble is negligibly small but did not provide an alternative explanation. It has been shown that during laser-induced plasma formation, water is partially dissociated into gaseous products. Atomic hydrogen and oxygen will largely recombine to form water but some molecular hydrogen and oxygen remain as long-lived gaseous products (Nikogosyan et al. Reference Nikogosyan, Oraevsky and Rupasov1983; Sato et al. Reference Sato, Tinguely, Oizumi and Farhat2013; Barmina et al. Reference Barmina, Simakin and Shafeev2016, Reference Barmina, Simakin and Shafeev2017).
Water dissociation can proceed through thermally driven hydrolysis (Mattsson & Desjarlais Reference Mattsson and Desjarlais2006, Reference Mattsson and Desjarlais2007) and by free-electron-mediated bond breaking. Liang, Zhang & Vogel (Reference Liang, Zhang and Vogel2019) showed that the average kinetic energy of free electrons (i.e. conduction band electrons) in luminescent plasmas in bulk water is 6.8 eV, and the high-energy tail of their energy spectrum reaches up to 14 eV. This energy is sufficient to break bonds by dissociative electron attachment (Cobut et al. Reference Cobut, Jay-Gerin, Frongillo and Patau1996; Fedor et al. Reference Fedor2006; Ram, Prabhudesai & Krishnakumar Reference Ram, Prabhudesai and Krishnakumar2009). When conduction band electrons solvate, they go through a process of ‘hydration’ until they are trapped at an energy level of 6.5 eV above the valence band (Linz et al. Reference Linz, Freidank, Liang, Vogelmann, Trickl and Vogel2015). The hydrated electrons may also contribute to dissociation and subsequent gas formation (Draganic & Draganic Reference Draganic and Draganic1971; Nikogosyan et al. Reference Nikogosyan, Oraevsky and Rupasov1983; Elles et al. Reference Elles, Shkrob, Crowell and Bradforth2007) because their energy (≥6.5 eV) is larger than the O–H bonding energy of § 5.2 eV (Toegel, Hilgenfeldt & Lohse Reference Toegel, Hilgenfeldt and Lohse2002; Maksyutenko, Rizzo & Boyarkin Reference Maksyutenko, Rizzo and Boyarkin2006). At lower plasma temperatures, free-electron-mediated gas generation dominates, while for T ≥ 3000 K thermal dissociation becomes the dominant mechanism (Lédé, Lapique & Villermaux Reference Lédé, Lapique and Villermaux1983; Mattsson & Desjarlais Reference Mattsson and Desjarlais2006). For T > 4500 K, most of the chemical bonds are broken, and only radicals of H and O are present (Lédé et al. Reference Lédé, Lapique and Villermaux1983; Jung, Jang & You Reference Jung, Jang and You2013). Correspondingly, the amount of non-condensable gas in laser-produced cavitation bubble was found to increase with growing plasma temperatures (Sato et al. Reference Sato, Tinguely, Oizumi and Farhat2013).
In the case of laser-induced bubble generation that is portrayed in detail in the present paper, the plasma temperature is relatively low (Tavg = 1550 K). Although it suffices to induce complete vaporization of the water within the plasma volume, it is so low that the bubble's non-condensable gas content must have been produced exclusively by free-electron-mediated water dissociation. Non-condensable gas can leave the bubble only by dissolution, which takes much longer than the lifetime of the transient cavitation bubble (Baffou et al. Reference Baffou, Polleux, Rigneault and Monneret2014). Therefore, the bubble size at late stages of its lifetime thus provides information on the amount of gas produced during breakdown. For the bubble of figure 6, we found a water vapour/gas ratio of approximately 2 : 1 during first collapse (§ 4.2.3), while in SBSL bubbles the vapour/gas ratio is 1 : 5 (Storey & Szeri Reference Storey and Szeri2000). In future, it will be interesting to investigate the influence of plasma energy density on water dissociation in more detail to elucidate its influence on the strength of the bubble collapse.
Most of the vapour produced during optical breakdown condenses during the initial adiabatic expansion but during the subsequent isothermal expansion phase up to Rmax, vapour again invades the bubble, and the vapour content of the maximally expanded bubble by far exceeds the amount of non-condensable gas (Storey & Szeri Reference Storey and Szeri2000). At collapse, vapour condenses again but during the final collapse stage, the bubble wall movement becomes so fast that some vapour is trapped inside the bubble. This process was modelled explicitly in several studies (Fujikawa & Akamatsu Reference Fujikawa and Akamatsu1980; Yasui Reference Yasui1995; Storey & Szeri Reference Storey and Szeri2000; Toegel et al. Reference Toegel, Gompf, Pecha and Lohse2000; Akhatov et al. Reference Akhatov, Lindau, Topolnikov, Mettin, Vakhitova and Lauterborn2001; Lauer et al. Reference Lauer, Hu, Hickel and Adams2012; Zein et al. Reference Zein, Hantke and Warnecke2013; Zhong et al. Reference Zhong, Eshragi, Vlachos, Dabiri and Ardekani2020). Unfortunately, the values of evaporation and condensation coefficients depend on pressure and temperature, and their values are still uncertain (Eames, Marr & Sabir Reference Eames, Marr and Sabir1997; Marek & Straub Reference Marek and Straub2001). We escaped the dilemma by using Rn as a fit parameter. In conjunction with (3.52)–(3.55), this approach is not just a makeshift but also an indirect way to determine the amount of vapour condensing during the bubble's expansion, collapse and rebound. Our results for the vapour mass reduction during first bubble collapse (table 1) are in good agreement with the findings of Akhatov et al. (Reference Akhatov, Lindau, Topolnikov, Mettin, Vakhitova and Lauterborn2001).
5.3. Extreme conditions during generation and collapse of laser-induced bubbles
Optical breakdown by tightly focused laser pulses produces plasmas with average volumetric energy densities reaching from a few kJ cm−3 to approximately 40 kJ cm−3, depending on laser pulse duration, wavelength, focusing angle and pulse energy (Vogel et al. Reference Vogel, Nahen, Theisen and Noack1996b). In this paper, we investigated the cavitation events produced by a weakly luminescent plasma with relatively small temperature and small non-condensable gas content that exhibits a particularly strong bubble collapse. The bubble collapse concentrates the potential energy of the expanded bubble at Rmax 1 into the volume of the collapsed bubble and the compressed liquid in its vicinity. For the vigorously collapsing bubble investigated in this paper, the volume ratio ${({R_{max1}}/{R_{min1}})^3}$ is 6.28 × 105. This energy concentration produces a peak pressure of 13.5 GPa and a temperature of 31 400 K. Under these conditions, luminescent plasma forms in the collapsed bubble (Baghdassarian et al. Reference Baghdassarian, Tabbert and Williams1999, Reference Baghdassarian, Chu, Tabbert and Williams2001; Ohl et al. Reference Ohl, Kurz, Geisler, Lindau and Lauterborn1999; Brenner et al. Reference Brenner, Hilgenfeldt and Lohse2002; Mattsson & Desjarlais Reference Mattsson and Desjarlais2006).
Upon collapse, most energy is carried away by shock wave emission and little energy remains for the rebounding bubble, which results in a large Rmax 1/Rmax 2 ratio. Figure 15 compares ratios from previous publications with values obtained in the present study. In previous studies with larger bubbles, buoyancy effects distorted the spherical shape, and the Rmax 1/Rmax 2 ratios were significantly smaller than in the present paper. Moreover, the relative importance of water vapour increases for larger bubbles because the amount of water contained in the expanded bubble scales with $R_{max}^3$, while the surface area through which condensing vapour can escape during collapse scales proportional to $R_{max}^2$ (Toegel et al. Reference Toegel, Gompf, Pecha and Lohse2000). The parameter combination explored in the present study (small bubble and relatively low plasma energy density) results in a particularly strong collapse.
How does the collapse of laser-induced bubbles compare with SBSL bubbles? With excitation in the kHz range, the most vigorous dynamics occurs around f ≈ 16 kHz (Toegel et al. Reference Toegel, Gompf, Pecha and Lohse2000) and with a driving pressure pa ≈ 0.145 MPa (Matula Reference Matula1999; Toegel & Lohse 2003). Simulations by Matula (Reference Matula1999) predicted a ratio Rmax/Rnc = 9.94 for a SBSL bubble driven at f ≈ 25 kHz and pa ≈ 0.142 MPa. Since a similar ratio (Rmax 1/Rnc 1 = 9.94) was found in the present study for a laser-induced bubble, similar collapse pressures are expected in both cases. Streak-photographic measurements of SBSL collapse pressure support this conclusion. Pecha & Gompf (Reference Pecha and Gompf2000) found a shock wave velocity of us = 4000 m s−1 at f ≈ 20 kHz and pa ≈ 0.139 MPa, and Weninger et al. (Reference Weninger, Evans and Putterman2000) reported us = 5930 m s−1 (Mach 4) for f ≈ 16.5 kHz and pa ≈ 0.145 MPa. An evaluation of these us data using the Hugoniot data of Rice & Walsh (Reference Rice and Walsh1957) yields ps = 5.3 GPa and ps = 15.4 GPa, respectively, close to the collapse pressure of 13.5 GPa obtained in the present paper.
In the context of SBSL, researchers pointed out that the wall of the collapsing bubble can launch an internal shock wave when its velocity exceeds the sound velocity inside the bubble (Roberts & Wu Reference Roberts and Wu1996; Lin & Szeri Reference Lin and Szeri2001; Brenner et al. Reference Brenner, Hilgenfeldt and Lohse2002). It was postulated that the geometrical focusing of the internal shock wave could produce a tiny spot in the bubble centre with strongly elevated pressure and temperatures up to 106 K (Wu & Roberts 1993). However, Storey & Szeri (Reference Storey and Szeri2000) demonstrated that endothermic reactions of water vapour trapped in the collapsing bubble significantly reduce peak temperature and pressure. An equilibrating factor for the pressure distribution within the bubble is the rapid increase of sound velocity upon collapse. The density of the bubble content can exceed the liquid density both for collapsing gas bubbles (Yuan et al. Reference Yuan, Ho, Chu and Leung2001) and vapour bubbles, and the sound velocity thus assumes much higher values than under ambient conditions. This effect is further enhanced by the rise of sound speed with increasing temperature (Vuong, Szeri & Young 1999). For a collapse pressure value of 13.5 GPa (figure 7), the corresponding sound velocity derived from the EOS data by Rice & Walsh (Reference Rice and Walsh1957) is 5400 m s−1. This value is much larger than the peak bubble wall velocity of 1793 m s−1, and the formation of an inner shock wave will thus be impeded. Lack of an inner shock wave justifies the assumption of a homogeneous bubble pressure made in the Gilmore model.
The pressure jump at the front of the external shock wave emitted upon breakdown and rebound is often so high that energy dissipation causes a temperature rise beyond the spinodal limit. This extends the region in which vaporization occurs after breakdown and upon rebound in a more effective way than heat conduction does (in table 2, ESW during the entire bubble life is 54.2 % of Eabs, and Econd amounts to 43 %). For the bubble's rebound phase, shock-wave-induced phase transitions have not yet been captured by time-resolved shadow or Schlieren photographs. However, in figure 5 it is visualized during the expansion of a high-density plasma produced by an energetic laser pulse. After plasma formation, reproducible timing of photographs is easier than for the rebound phase immediately after collapse, and the size of the affected region is large enough to be resolved by optical imaging. During the rebound of a collapsed spherical bubble, similar processes will occur. In figure 10(b), a shock front exhibiting a maximum pressure jump of ≈8 GPa develops, which results in a temperature jump to 436 °C (Rice & Walsh Reference Rice and Walsh1957) that produces a phase transition and lowers surface tension and viscous damping. While models for heat and mass transfer by evaporation and condensation at the bubble wall and heat conduction are already available, the ‘convective’ heat transport by the shock wave has not yet been considered in any bubble model. Its inclusion remains a challenge for the future.
5.4. Pressure decay during shock wave propagation
The first studies on acoustic emission during spherical bubble collapse by Hickling & Plesset (Reference Hickling and Plesset1964) and Akulichev et al. (Reference Akulichev, Boguslav, Ioffe and Naugolny1968) as well as the simulations by Fujikawa & Akamatsu (Reference Fujikawa and Akamatsu1980) were based on assumptions about the bubble's gas content rather than on measured R(t) curves or oscillation times. Akhatov et al. (Reference Akhatov, Lindau, Topolnikov, Mettin, Vakhitova and Lauterborn2001) presented p(r) curves for a model fit to experimental R(t) data but undertook no systematic study of the shock wave's pressure decay. The present paper presents the first evidence-based simulation of shock wave formation and pressure decay after a vigorous collapse of spherical laser-induced cavitation bubbles.
Our results differ strongly from earlier results for lower collapse pressure, which yielded a pressure decay resembling that of acoustic transients with ${p_{peak}} \propto {r^{ - 1}}$ (Hickling & Plesset Reference Hickling and Plesset1964; Fujikawa & Akamatsu Reference Fujikawa and Akamatsu1980; Koch et al. Reference Koch, Lechner, Reuter, Kohler, Mettin and Lauterborn2016). We observe a rapid formation of shock waves as previously reported by Akulichev et al. (Reference Akulichev, Boguslav, Ioffe and Naugolny1968) and Ebeling (Reference Ebeling1978), together with a pressure decay as fast as ${p_{peak}} \propto {r^{ - 1.75}}$. A similar decay rate, with a maximum slope in the log–log plot of −1.79, has been obtained in previous simulations for the breakdown shock wave produced by a 10 mJ ns laser pulse that was emitted from high-density plasma with 8.8 GPa initial pressure (Vogel et al. Reference Vogel, Busch and Parlitz1996a).
Experimentally observed pressure decays are usually steeper than predictions by the Gilmore model, with steepest slopes ≥2.0 (Vogel et al. Reference Vogel, Busch and Parlitz1996a; Lai et al. Reference Lai, Geng, Zheng, Yao, Zhong and Wang2021). This is because the evaluation of experimental data is based on Hugoniot data valid up to 25 GPa, while the Gilmore model employs the Tait EOS that for very large shock wave velocities yields pressure values that are too low (see § 5.19).
It is interesting to note that the experimentally observed pressure decay measured by Vogel et al. (Reference Vogel, Busch and Parlitz1996a) in the direction perpendicular to the laser beam axis was initially weaker and then faster than the decay predicted by the simulations for spherical bubble dynamics. This discrepancy is likely caused by the elongated plasma shape in the experiments, which introduces a cylindrical component in the near-field shock wave that goes along with a slower pressure decay in the direction perpendicular to the laser beam axis (Schoeffmann, Schmidt Kloiber & Reichel Reference Schoeffmann, Schmidt Kloiber and Reichel1988). However, the shock wave propagation in the far field exhibits radial symmetry even when the breakdown region is elongated (Tagawa et al. Reference Tagawa, Yamamoto, Hayasaka and Kameda2016). Therefore, the pressure decay must be faster in the transition zone between near and far field to make up for the slower near-field decay. A similar phenomenon is seen in figure 13 for the velocity of bubble expansion around elongated plasmas. Future experiments with laser-induced bubbles exhibiting minimum deviations from spherical shape will enable a more precise comparison with numerical simulations than possible to date.
5.5. Energy partitioning
In this paper, we established a complete energy balance for laser-induced spherical cavitation bubbles, while previous studies focused on the energy partitioning in bubble and shock wave energy (Vogel et al. Reference Vogel, Nahen, Theisen and Noack1999b; Tinguely et al. Reference Tinguely, Obreschkow, Kobel, Dorsaz, de Bosset and Farhat2012). The new approach cannot only quantify the total amount of shock wave energy emitted after breakdown and collapse but also distinguish between fractions of the collapse shock wave originating from energy stored in the compressed bubble and from the compressed liquid surrounding it. This way we revealed strong differences between the energy partitioning after optical breakdown and bubble collapse. During breakdown, the laser energy is stored inside the plasma, and a compression wave affects the surrounding liquid only after the plasma has started to expand (figure 9). Under these conditions, a relatively large fraction of the plasma energy can be converted into bubble energy. With moderate plasma energy density, 59.2 % of the absorbed energy was transformed into mechanical energy $(E_{SW}^{bd} + E_{pot}^{max1} + {W_{visc}} + U_{int}^{max1})$, and from this fraction 55.4 % went into bubble energy and 38.9 % into shock wave emission (table 2). By contrast, during the rebound after the first bubble collapse, less than 2 % was transformed into bubble energy, and 96.4 % of the energy stored in the compressed bubble content and liquid was radiated away acoustically. This is because both the bubble content and the surrounding liquid are compressed during collapse (see figure 10), with the energy content of the compressed liquid being much larger than that of the bubble (table 2). Therefore, most energy is radiated away acoustically upon rebound and only a small fraction originating from the internal energy of the compressed bubble content can contribute to bubble formation. For the bubble investigated in detail in the present paper, the ratio of the energy contributions from compressed liquid and bubble content was ESWL/ESWB = 18.45.
These findings agree qualitatively with the picture on acoustic emission after bubble collapse obtained with models based on solutions of the Navier–Stokes equations (Fuster et al. Reference Fuster, Dopazo and Hauke2011). A quantitative comparison between their results and the present results is difficult because the collapse pressures differ significantly. In our case, the collapse pressure was 13.5 GPa, whereas the largest collapse pressure investigated by Fuster et al. (Reference Fuster, Dopazo and Hauke2011) was only 1.5 GPa according to the Gilmore model, and 2.6 GPa according to their full model. The ratio ESWL/ESWB increases with decreasing gas content of the collapsing bubble, when less internal energy is stored in the bubble itself and more in the surrounding liquid. This goes along with increasing collapse pressure and stronger acoustic emission from the liquid surrounding the bubble. Therefore, ESWL/ESWB is particularly high in the present paper.
With increasing plasma energy density, a smaller fraction of the absorbed laser energy is required for vaporization of the plasma volume and an ever-larger percentage is converted into mechanical energy (Vogel et al. Reference Vogel, Nahen, Theisen and Noack1999b). Furthermore, an ever-larger part of the mechanical energy appears as shock wave energy (Lai et al. Reference Lai, Geng, Zheng, Yao, Zhong and Wang2021). For a 6 ns, 10 mJ pulse with ε = 40 kJ cm−3, the shock wave energy was found to be more than two times larger than the bubble energy (Vogel et al. Reference Vogel, Nahen, Theisen and Noack1999b), different from the case of figure 6 in the present paper, where for ε = 8.7 kJ cm−3 the shock wave energy is slightly smaller than the bubble energy. This finding warrants future detailed investigations of the dependence of energy partitioning on plasma energy density.
Figure 15 shows that for very small bubble sizes, the Rmax 1/Rmax 2 ratio increases strongly with decreasing bubble size, i.e. more energy is dissipated upon bubble collapse. This change is most likely related to the increasing role of surface tension and viscosity with decreasing bubble size as indicated by the 1/R proportionality of the last two terms of (3.3). The increase of the Laplace pressure arising from surface tension enhances the vigour of the collapse while, at the same time, an ever-larger part of the deposited energy is dissipated by viscous damping. The tools presented in this paper enable a quantitative investigation of the changes in bubble dynamics and energy partitioning for ${R_{max}} \to 0$.
6. Conclusions
We established a hybrid experimental/simulation approach for providing a rapid and comprehensive characterization of laser-induced bubble oscillations and shock wave emission. The experimental part consists of a photographic characterization of the size of the laser-induced plasma and a single-shot probe beam scattering method for recording the bubble oscillation times with high temporal resolution. The excellent time resolution, large dynamic range and high sensitivity of the method enables us to cover the entire bubble lifetime from early large-amplitude nonlinear oscillations to late oscillations of the residual gas bubble with sub-nanometre amplitude.
Simulations are performed based on the Gilmore model with a van der Waals hard core that has been extended by a description of laser-induced bubble formation considering the shock-wave-induced jump start of the bubble wall, an automated determination of the shock front location for pressure transients with very large amplitude, and an energy balance encompassing the entire bubble lifetime. The results of the experiments and calculations complement each other and yield a rich and detailed picture of the events during spherical laser-induced cavitation bubble oscillations.
Laser-induced bubble dynamics is a microlaboratory for high-pressure/density/temperature hydrodynamics, plasma physics and chemistry that enables us to study a large number of nonlinear phenomena and extreme states in a tabletop environment. Laser-induced bubbles offer the option to study both spherical and aspherical bubble dynamics under controlled conditions and are of great practical importance in biophotonics and biomedicine as well as in laser ablation in liquids. Their investigation will continue to provide fruitful insights if modelling and experimental tools are further advanced. However, the challenges for experimental coverage with high spatial and temporal resolution are extremely high because the collapse times show much larger shot-to shot fluctuations than the oscillation times in SBSL.
Modelling is also very demanding because of the simultaneous occurrence of a multitude of nonlinear phenomena and additional challenges posed by aspherical dynamics. Volume of fluid methods are a versatile tool providing new insights, especially for aspherical dynamics. However, for a better understanding of spherical bubble dynamics in the context of biomedical applications, the main focus lies on the dependence of the dynamics on laser parameters and properties of the breakdown medium. Our spherical bubble model combines relative simplicity with large information content based on few readily available experimental data, which makes it a useful tool for evidence-based investigations of parameter dependencies. In this paper, we demonstrated the large potential of this approach on one example, where the bubble dynamics could be traced through more than 100 oscillations. In future work, this tool will be applied to characterize the changes in bubble dynamics and energy partitioning in dependence on plasma energy density and for ${R_{max}} \to 0$.
Supplemental material
Supplementary materials are available at https://doi.org/10.1017/jfm.2022.202.
Funding
This work was supported by U.S. Air Force Office of Scientific Research (X.-X.L, A.V, grant numbers FA9550-15-1-0326, FA9550-18-1-0521).
Declaration of interest
The authors report no conflict of interest.