Comprehensive analysis of spherical bubble oscillations and shock wave emission in laser-induced cavitation

Abstract

The dynamics of spherical laser-induced cavitation bubbles in water is investigated by plasma photography, time-resolved shadowgraphs and sensitive single-shot probe beam scattering that portrays the transition from initial nonlinear to late linear oscillations. The frequency of late oscillations yields the bubble's gas content. Numerical simulations with an extended Gilmore model using plasma size as input and oscillation times as fit parameter provide insights into experimentally not accessible bubble parameters and shock wave emission. Model extensions include a term covering the initial shock-driven acceleration of the bubble wall, an automated method determining shock front position and pressure decay and a complete energy balance for the partitioning of absorbed laser energy into vaporization, bubble and shock wave energy and dissipation through viscosity and condensation. These tools are used for analysing a scattering signal covering 102 oscillation cycles from a bubble with 36 μm maximum radius produced by a plasma with 1550 K average temperature. Predicted bubble wall velocities during expansion agree well with experimental data. Upon first collapse, most energy was stored in the compressed liquid around the bubble and radiated away acoustically. The collapsed bubble contained more vapour than gas and had a pressure of 13.5 GPa. The decay of the rebound shock wave pressure with radius r was initially $\mathrm{\ \propto }{r^{ - 1.8}}$, and energy dissipation at the shock front heated the liquid near the bubble wall to temperatures above the superheat limit. The shock-induced temperature rise reduces damping during late bubble oscillations. Damping in first collapse increases significantly for small bubbles with less than 10 μm radius.

1. Introduction

Laser-induced cavitation plays an important role in many fields, especially in laser materials processing in liquid environments (Barcikowski et al. 2019; Kanitz et al. 2019) and in biomedicine and biophotonics, where it enables precise surgery on cells and within transparent tissues (Vogel et al. 1990; Juhasz et al. 1999; Koenig et al. 1999; Tirlapur & Koenig 2002; Yanik et al. 2004; Vogel et al. 2005; Chung & Mazur 2009; Palanker et al. 2010; Stevenson et al. 2010; Hoy et al. 2014). 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. 1996b; Lauterborn & Vogel 2013). 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 1966; Lauterborn & Bolle 1975; Tomita & Shima 1986; Vogel, Lauterborn & Timm 1989; Vogel, Busch & Parlitz 1996a; Philipp & Lauterborn 1998; Baghdassarian, Tabbert & Williams 1999; Brujan et al. 2001; Lindau & Lauterborn 2003; Brujan & Vogel 2006; Obreschkow et al. 2011; Obreschkow et al. 2013; Reuter et al. 2017; Sagar & el Moctar 2020; Podbevsek et al. 2021). 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 1977; Lauterborn & Kurz 2010). 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 2010). 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 1952) 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 ₀ for the simulations, and the duration of the subsequent bubble oscillations, T_osci, 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 1974; Müller et al. 2009; Kroeninger et al. 2010; Obreschkow et al. 2013). 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. 1996a; Schaffer et al. 2002; Venugopalan et al. 2002). 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. 2008). 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. 1997; Matula 1999; Gompf & Pecha 2000; Weninger, Evans & Putterman 2000; Schaffer et al. 2002). In the present paper, we use the confocal scheme introduced by Vogel et al. (2008) 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. 2002; Vogel et al. 2008; Obreschkow et al. 2013; Sinibaldi et al. 2019). Additionally, buoyancy effects that could lead to a movement out of the probe beam focus and to jet formation upon collapse (Benjamin & Ellis 1966; Zhang et al. 2015) 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 1966; Blake & Gibson 1987) or under zero gravity conditions during parabolic flights (Obreschkow et al. 2011, 2013). 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 (1980) 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 1999). 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⁻¹ are observed (Vogel et al. 1996a). 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. 1996b) and the contribution of the particle velocity behind the detaching shock front producing a jump start of the bubble wall velocity. Gilmore (1952) already presented a first-order approximation for the quick start of the wall movement of a bubble starting to collapse from R = R_max 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 1980; Akhatov et al. 2001; Lauer et al. 2012; Zein, Hantke & Warnecke 2013; Peng et al. 2020; Zhong et al. 2020; Aganin & Mustafin 2021) 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 2001), making predictions uncertain. Akhatov et al. (2001) 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, R_nbd, and during collapse, R_{nc 1} and R_{nc 2}, as fit parameters, which are chosen such that model predictions match the observed oscillation times for the first three oscillations, T_{osc 1} to T_{osc 3}, following Vogel et al. (1996a), Lauterborn & Kurz (2010) and Koch et al. (2016). 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 R_{max 1} is obtained by assuming equilibrium vapour pressure at room temperature, and the size the residual bubble containing mainly non-condensable gas, R_res, is derived from the frequency of the late linear oscillations. Comparison of R_{nc 1} and R_res 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 ₀ and T_osci 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. 1999a; Vogel et al. 2008; Obreschkow et al. 2011, 2013; Sinibaldi et al. 2019). 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 1971; Prosperetti & Hao 1999; Yuan et al. 2001). 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 1966). 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

(2.1)\begin{equation}\delta = \sqrt {\frac{{{\rho _0}g{R_{max}}}}{{{p_\infty } - {p_v}}}} ,\end{equation}

which expresses the ratio of bubble collapse time to the time it takes an inviscid bubble of radius R_max to rise one radius from rest driven by buoyancy forces (Blake, Taib & Doherty 1986; Best & Kucera 1992). Here, ${p_\infty }$ denotes the ambient pressure, p_v the vapour pressure at ambient conditions, ρ ₀ 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 R_max = 100 μm (δ = 0.003), in accordance with previous observations of stable sonoluminescence for bubbles of ≈90 μm maximum radius (Gompf & Pecha 2000; Weninger et al. 2000). Investigation of small bubbles not only minimizes the influence of buoyancy but also provides a large dimensionless stand-off distance γ = d/R_max from the front lens of the microscope objective. For R_max = 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 R_max/T_osc (Kroeninger et al. 2010; Obreschkow et al. 2013). Therefore, it is sufficient to measure T_osc to characterize spherical bubble oscillations.

Figure 1. Experimental arrangements for the investigation of the behaviour of spherical laser-induced bubbles in water. (a) Set-up with confocal adjustment of three microscope objectives enabling us to generate highly spherical bubbles, record their oscillations with a cw probe laser beam and take high-resolution images of plasma luminescence. (b) Illustration of the directly transmitted and multiply reflected parts of the probe laser beam that interfere behind large bubbles. Light scattering and interference are detected by the AC-coupled photoreceiver in (a) and recorded using a digital oscilloscope. (c) Set-up for time-resolved photography of bubble formation and shock wave emission at times up to t = 120 ns. An optically delayed frequency-doubled portion of the pump laser pulse is used for illumination. Hydrophone signals of breakdown and collapse shock waves are recorded to monitor the bubble oscillation time for each shot.

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. (2016). 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 T_tra 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. 2008; Linz et al. 2016). 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, T_osc (Vogel et al. 2008). 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 R_max = 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 < R_max < 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 1952; Lauterborn & Kurz 2010) 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 1948). 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 1993; Lauterborn & Kurz 2010).

The bubble dynamics is described by the equation

(3.1)\begin{equation}\left( {1 - \frac{U}{C}} \right)R\dot{U} + \frac{3}{2}\left( {1 - \frac{U}{{3C}}} \right){U^2} = \left( {1 + \frac{U}{C}} \right)H + \frac{U}{C}\left( {1 - \frac{U}{C}} \right)R\frac{{\textrm{d}H}}{{\textrm{d}R}}.\end{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

(3.2)\begin{equation}H = \int_{p{|_{r \to \infty }}}^{p{|_{r = R}}} {\frac{{\textrm{d}p{\kern 1pt} (\rho )}}{\rho }} ,\end{equation}

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

(3.3)\begin{equation}P = p\left|{_{r = R} = \left( {{p_\infty } + \frac{{2\sigma }}{{{R_n}}}} \right)} \right.{\left( {\frac{{R{\kern 1pt}_n^3 - R{\kern 1pt}_{v\,\textrm{d}W}^3}}{{{R^3} - R{\kern 1pt}_{v\,\textrm{d}W}^3}}} \right)^\kappa } - \frac{{2\sigma }}{R} - \frac{{4\mu }}{R}U,\end{equation}

where σ denotes the surface tension, μ the dynamic shear viscosity and κ the ratio of the specific heat at constant pressure and volume. The symbol R_n 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 1988)

(3.4)\begin{equation}\frac{{p + B}}{{{p_\infty } + B}} = {\left( {\frac{\rho }{{{\rho_\infty }}}} \right)^n},\end{equation}

which leads to the following relationships for the sound velocity C and enthalpy H at the bubble wall:

(3.5)\begin{gather}C = \sqrt {c_\infty ^2 + (n - 1)H} ,\end{gather}

(3.6)\begin{gather}H = \frac{{n({p_\infty } + B)}}{{(n - 1){\rho _\infty }}}\left[ {{{\left( {\frac{{P + B}}{{{p_\infty } + B}}} \right)}^{(n - 1)/n}} - 1} \right],\end{gather}

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.7)\begin{equation}\frac{{\textrm{d}H}}{{\textrm{d}R}} = \frac{1}{{{\rho _0}}}{\left( {\frac{{{p_\infty } + B}}{{P + B}}} \right)^{1/n}} \times \left( { - 3\kappa {R^2}\left( {{p_\infty } + \frac{{2\sigma }}{{{R_n}}}} \right)\frac{{{{(R_n^3 - R_{v\,\textrm{d}W}^3)}^\kappa }}}{{{{({R^3} - R_{v\,\textrm{d}W}^3)}^{\kappa + 1}}}} + \frac{{2\sigma }}{{{R^2}}} + \frac{{4\mu U}}{{{R^2}}}} \right).\end{equation}

3.2. Description of laser-induced bubble initiation

Following Vogel et al. (1996a), 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 ₀, 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 R_n from its small initial value R_n = R ₀ at the beginning of the pulse to a much larger final value R_nbd. 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 R_n 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 R_n value beyond R ₀ implies that the pressure inside the bubble rises and that the bubble starts to expand. Iteratively, we determine the R_nbd value for which the calculation yields the same oscillation time T_{osc 1} or maximum radius R_{max 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 P_L during the pulse is modelled by a sin² function with duration τ_L (full-width at half-maximum) and total duration 2τ_L

(3.8)\begin{equation}{P_L} = {P_{L0}}{\sin ^2}\left( {\frac{\rm \pi}{{2{\tau_L}}}t} \right),\quad 0 \le t \le 2{\tau _L}.\end{equation}

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 E_L absorbed up to this time, Vogel et al. (1996a) derived an equation for the temporal development of the equilibrium radius R_n during the laser pulse

(3.9)\begin{equation}{R_n}(t) = {\left\{ {R_0^3 + \frac{{R_{nbd}^3 - R_0^3}}{{2{\tau_L}}}\left[ {t - \frac{{{\tau_L}}}{\rm \pi}\sin \left( {\frac{\rm \pi}{{{\tau_L}}}t} \right)} \right]} \right\}^{1/3}}.\end{equation}

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 ₀ equals the initial particle velocity u_p behind the shock front, whereby the shock pressure is identical with the initial plasma and bubble pressure, p_s = P.

Gilmore considered the case where the internal bubble pressure, P_i, 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

(3.10)\begin{equation}{U_0} = \int_0^H {\frac{{\textrm{d}h}}{c}} \approx \frac{H}{C} \approx \frac{{P - {p_\infty }}}{{{\rho _\infty }{c_\infty }}}.\end{equation}

The approximate expression is accurate when $|H|\ll {C^2}$, which for water corresponds to $|P_{i} - {p_\infty }|\ll 2000\ \textrm{MPa}$ (Gilmore 1952). 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 ₀ based on the Hugoniot curve data from Rice & Walsh (1957). 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

(3.11)\begin{equation}{u_p} = {c_1}({10^{({u_s} - {c_\infty })/{c_2}}} - 1),\end{equation}

where u_s is the shock wave velocity and the constants are c ₁ = 5190 m s⁻¹, c ₂ = 25 306 m s⁻¹ and c ₀ is the sound velocity, c_∞ = 1483 m s⁻¹. By rearranging (3.11), u_s can be expressed as a function of u_p

(3.12)\begin{equation}{u_s} = {c_\infty } + {c_2}{\log _{10}}({u_p}/{c_1} + 1).\end{equation}

Using (3.11) and the conservation of momentum at a shock front, ${p_s} - {p_\infty } = {u_s}{u_p}{\rho _\infty }$ (Duvall & Fowles 1963), one can link p_s to u_s

(3.13)\begin{equation}{p_s} = {c_1}{\rho _\infty }{u_s}({10^{({u_s} - {c_\infty })/{c_2}}} - 1) + {p_\infty },\end{equation}

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 p_s and u_p

(3.14)\begin{equation}{p_s} = {\rho _\infty }{u_p}[{c_\infty } + {c_2}{\log _{10}}({u_p}/{c_1} + 1)] + {p_\infty }.\end{equation}

If ${u_p} \ll {c_1}$, the second term in the bracket can be dropped, which leads to

(3.15)\begin{equation}{p_s} = {\rho _\infty }{u_p}{c_\infty } + {p_\infty }.\end{equation}

Immediately after breakdown, u_p = U ₀ and p_s = P. After resolving (3.15) for U ₀, we get

(3.16)\begin{equation}{u_p} = {U_0} = (P - {p_\infty })/{\rho _\infty }{c_\infty },\end{equation}

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

(3.17)\begin{equation}\dot{U} ={-} \frac{3}{2}\frac{{{U^2}}}{R}\frac{{C - U/3}}{{C - U}} + \frac{H}{R}\frac{{C + U}}{{C - U}} + \frac{U}{C}\frac{{\textrm{d}H}}{{\textrm{d}R}} + {\dot{u}_p}.\end{equation}

The term ${\dot{u}_p}$ is derived from (3.16) as

(3.18)\begin{equation}{\dot{u}_p} = \left\{ {\begin{array}{*{20}{@{}l}} {\dfrac{{\dot{P}{\kern 1pt} }}{{{\rho_\infty }{c_\infty }}}}&{\textrm{for}\;0 \le t \le 2{\tau_L},}\\ 0&{\textrm{otherwise}\textrm{.}} \end{array}} \right.\end{equation}

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 ₀. 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

(3.19)\begin{equation}P = \frac{{{p_\infty }}}{{R_0^4}}R_n^4(t) + \frac{{2\sigma }}{{R_0^4}}R_n^3(t),\end{equation}

with R_n(t) given by (3.9). The time derivative of (3.19) reads

(3.20)\begin{equation}\dot{P} = \frac{{4{p_\infty }{R_n}(t) + 6\sigma }}{{R_0^4}}R_n^2(t){\dot{R}_n}(t),\end{equation}

and the time derivative of (3.9) is

(3.21)\begin{equation}{\dot{R}_n}(t) = \frac{1}{3}R_n^{ - 2}(t)\frac{{R_{nbd}^3 - R_0^3}}{{2{\tau _L}}}\left[ {1 - \cos \left( {\frac{\rm \pi}{{{\tau_L}}}t} \right)} \right].\end{equation}

By inserting (3.21) into (3.20) one gets

(3.22)\begin{equation}\dot{P} = \frac{{2{p_\infty }{R_n}(t) + 3\sigma }}{{3{\tau _L}R_0^4}}(R_{nbd}^3 - R_0^3)\left[ {1 - \cos \left( {\frac{\rm \pi}{{{\tau_L}}}t} \right)} \right],\end{equation}

and by inserting (3.22) into (3.18) one finally obtains

(3.23)\begin{equation}{\dot{u}_p} = \frac{{2{p_\infty }{R_n}(t) + 3\sigma }}{{3R_0^4{\tau _L}{\rho _\infty }{c_\infty }}}(R_{nbd}^3 - R_0^3)\left[ {1 - \cos \left( {\frac{\rm \pi}{{{\tau_L}}}t} \right)} \right],\end{equation}

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 P_max = 1.31 GPa, and the corresponding particle velocity is only approximately $500\ \textrm{m}\;{\textrm{s}^{ - 1}}$ (Rice & Walsh 1957).

Figure 2. Time evolution of bubble wall velocity during the early expansion phase for three modelling approaches: (i) particle velocity behind the shock wave front is not considered (no jump start); (ii) the first-order approximation of ${\dot{u}_p}$ in (3.23) is used to consider the evolution of particle velocity during the laser pulse; (iii) the second-order approximation in (3.30) is used. For all simulations, the input parameters are R ₀ = 1.33 μm and R_nbd = 13.88 μm, which correspond to the signal in figure 9 that will later be analysed in detail.

Besides overestimating U ₀, 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 u_p in (3.14) and formulate a second-order approximation considering the second term in the bracket through its Taylor expansion

(3.24)\begin{equation}{\log _{10}}\left( {\frac{{{u_p}}}{{{c_1}}} + 1} \right) = \frac{{{u_p}}}{{\log (10){c_1}}} - \frac{{u_p^2}}{{\log (100)c_1^2}} + \frac{{u_p^3}}{{\log (1000)c_1^3}} - \cdots .\end{equation}

If u_p 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

(3.25)\begin{equation}P = {\rho _\infty }{c_\infty }{u_p} + \frac{{{\rho _\infty }{c_2}}}{{\log (10){c_1}}}u_p^2 + {p_\infty }.\end{equation}

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)$

(3.26)\begin{equation}\underbrace{{\frac{{{\rho _\infty }{c_2}}}{{\log (10){c_1}}}}}_{A}u_p^2 + \underbrace{{{\rho _\infty }{c_\infty }}}_{B}{u_p} - P = 0.\end{equation}

For A, B, P > 0, such equations have a positive real and a negative imaginary root

(3.27)\begin{equation}{u_p} = \frac{{ - B \pm \sqrt {{B^2} + 4AP} }}{{2A}}.\end{equation}

Inserting A and B in (3.26) into the positive root of (3.27), we obtain

(3.28)\begin{equation}{u_p} = \frac{{\sqrt {\rho _\infty ^2c_\infty ^2 + 4\frac{{{\rho _\infty }{c_2}}}{{\log (10){c_1}}}P} - {\rho _\infty }{c_\infty }}}{{\frac{{2{\rho _\infty }{c_2}}}{{\log (10){c_1}}}}}.\end{equation}

The time derivative of this equation is

(3.29)\begin{equation}{\dot{u}_p} = \frac{{\dot{P}}}{{\sqrt {\rho _\infty ^2c_\infty ^2 + \frac{{4{\rho _\infty }{c_2}}}{{\log (10){c_1}}}P} }}.\end{equation}

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

(3.30)\begin{equation}{\dot{u}_p} = \frac{1}{{\sqrt {\rho _\infty ^2c_\infty ^2 + \frac{{4{\rho _\infty }{c_2}}}{{\log (10){c_1}}}P} }}\frac{{2{p_\infty }{R_n}(t) + 3\sigma }}{{3R_0^4{\tau _L}}}(R_{nbd}^3 - R_0^3)\left[ {1 - \cos \left( {\frac{\rm \pi}{{{\tau_L}}}t} \right)} \right],\end{equation}

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 R_n(t) from (3.9) was used to calculate the pressure distribution in the liquid surrounding the cavitation bubble (Gilmore 1952; Knapp, Daily & Hammitt 1970). 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 1948). The Kirkwood–Bethe hypothesis leads to the differential equations

(3.31)\begin{gather}\dot{u} ={-} \frac{1}{{c - u}}\left[ {(c + u)\frac{y}{{{r^2}}} - \frac{{2{c^2}u}}{r}} \right],\quad \dot{r} = u + c,\end{gather}

(3.32)\begin{gather}\textrm{with}\;c = {c_\infty }{\left( {\frac{{p + B}}{{{p_\infty } + B}}} \right)^{(n - 1)/2n}}.\end{gather}

The pressure p at $r = r(t)$ is given by

(3.33)\begin{equation}p = ({p_\infty } + B){\left[ {\left( {\frac{y}{r} - \frac{{{u^2}}}{2}} \right) \cdot \frac{{(n - 1){\rho_\infty }}}{{n({p_\infty } + B)}} + 1} \right]^{n/(n - 1)}} - B.\end{equation}

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 1977; Landau & Lifschitz 1987). 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 1947; Cole 1948; Rice & Walsh 1957; Duvall & Fowles 1963; Müller 2007).

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 R_nc is the equilibrium radius of the bubble relevant for the collapse phase. It is considerably smaller than R_nbd immediately after optical breakdown because most of the water vapour produced during bubble generation condenses during the oscillation (Ebeling 1978). The R_nc value is chosen such that the calculation yields the same oscillation time of the rebounding bubble, T_{osc 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 2006, 2007), resulting in diatomic molecules such as H₂ and O₂, which have a larger adiabatic exponent of 1.4 at room temperature (Fujikawa & Akamatsu 1980). 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 1999; Brenner, Hilgenfeldt & Lohse 2002). Therefore, some researchers use a continuously changing time-dependent value of κ (Brenner et al. 2002), and others switch from the isothermal to the adiabatic value, when the bubble radius passes the equilibrium radius (Barber et al. 1997; Yuan et al. 2001). Nevertheless, following Lauterborn & Kurz (2010), 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. 1996; Mattsson & Desjarlais 2006; Elles et al. 2007; Müller et al. 2009). 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 1983; Barmina, Simakin & Shafeev 2016, 2017). 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. 2002), diffusion of dissolved gas into a laser-induced cavitation bubble is negligibly small (Akhatov et al. 2001). The degree of water dissociation depends on temperature and, thus, on the initial plasma energy density (Mattsson & Desjarlais 2006; Sato et al. 2013). 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, R_res. The use of different R_n 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 1978). The R_nbd, R_{nc 1} and R_{nc 2} values are chosen by fitting the predicted bubble dynamics to measured values of T_{osc 1}, T_{osc 2} and T_{osc 3}, respectively. After the second collapse, the R_n 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 R_n 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 = R_max 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 1980; Prosperetti & Hao 1999; Akhatov et al. 2001). 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 R_n 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 R_n is gradual and not stepwise as assumed in our simulations, where R_n is reduced at R_{max 1} and R_{max 2}, but the stepwise reduction does not influence the predicted dynamics.

Fujikawa & Akamatsu (1980) and Yasui (1995) 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 R_nc 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 = R_max starting at room temperature with a virtual bubble pressure ${p_{R\,max,virt}}$ corresponding to the amount of gas represented by R_nc. This virtual starting pressure is lower than the real pressure given by the equilibrium vapour pressure at R = R_max. For determining ${p_{R\,max,virt}}$, we must consider that R_nc 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

(3.34)\begin{equation}{p_{R\,max,virt}} = {p_\infty }{\left( {\frac{{{R_{nc}}}}{{{R_{max}}}}} \right)^3}.\end{equation}

For an adiabatic collapse, pressure and temperature are linked by

(3.35)\begin{equation}p_1^{(1 - \kappa )}T_1^\kappa = p_2^{(1 - \kappa )}T_2^\kappa ,\end{equation}

which provides

(3.36)\begin{equation}{T_{coll}} = 293\ \textrm{K}{\left( {\frac{{{p_{R\,max,virt}}}}{{{p_{coll}}}}} \right)^{(1 - \kappa )/\kappa }},\end{equation}

for the collapse temperature. With κ = 4/3 and by inserting (3.34) into (3.36) we get

(3.37)\begin{equation}{T_{coll}} = 293\ \textrm{K}{\left( {\frac{{{p_{coll}}}}{{{p_\infty }}}} \right)^{1/4}}{\left( {\frac{{{R_{max}}}}{{{R_{nc}}}}} \right)^{3/4}}.\end{equation}

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. (2001) 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 2010)

(3.38)\begin{equation}{\nu _0} = \frac{1}{{2{\rm \pi}{R_{nres}}\sqrt {{\rho _\infty }} }}\sqrt {3\kappa \left( {{p_\infty } + \frac{{2\sigma }}{{{R_{nres}}}} - {p_v}} \right) - \frac{{2\sigma }}{{{R_{nres}}}} - \frac{{4{\mu ^2}}}{{{\rho _\infty }R_{nres}^2}}} .\end{equation}

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 R_nc 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 1948; Ebeling 1978). 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 1988; Vogel et al. 1996a; Vogel et al. 1999b) focused on an experimental investigation of energy partitioning for millimetre-sized bubbles by measuring the bubble's potential energy

(3.39)\begin{equation}{E_{pot}} = (4/3){\rm \pi} R_{max}^3({p_\infty } - {p_v}),\end{equation}

and shock wave energy

(3.40)\begin{equation}{E_{SW}} = \frac{{4{\rm \pi} R_m^2}}{{{\rho _\infty }{c_\infty }}}\int {p_s^2\,\textrm{d}t} ,\end{equation}

with R_m denoting the distance of the measurement location from the emission centre. Measurement of the temporal shock wave profile needed for the determination of E_SW is challenging (Vogel & Lauterborn 1988; Tinguely et al. 2012; Lauterborn & Vogel 2013), 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 1948; Vogel et al. 1996a; Vogel et al. 1999b). 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. 1996a). Another alternative was to determine the energy dissipation from the decay of shock wave pressure with propagation distance that was obtained by measuring u_s(r) (Vogel et al. 1999b). 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. (2012) 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 R_{max 1} and R_{max 2}). Considering the change of internal energy ΔU_int 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 ΔU_int 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, E_abs. During optical breakdown, the energy absorbed in the plasma volume ${V_P} = (4/3){\rm \pi} R_0^3$ partitions into vaporization energy

(3.41)\begin{equation}{E_v} = {\rho _\infty }(4/3){\rm \pi} R_0^3[{C_p}({T_2} - {T_1}) + {L_V}],\end{equation}

and an internal energy gain ΔU_int of the heated, pressurized gas volume. Here, T ₁ and T ₂ denote the room temperature (20°) and boiling temperature of water (100 °C), respectively, C_p = 4187 J (K kg)⁻¹ is the isobaric heat capacity of water at 20 °C and L_V = 2256 kJ kg⁻¹ is the latent heat of vaporization at 100 °C. An equation for ΔU_int will be given further below.

Figure 3. Energy partitioning pathways for laser-induced cavitation bubbles separated into four phases: bubble expansion, first collapse, first rebound and afterbounces. The respective phases are indicated by the superscripts $exp$, $coll$, $reb$ and $res$ in the symbols denoting the energy fractions. Solid arrows indicate the conversion of energy fractions that enter the next bubble oscillation phase, which include the vaporization energy ${E_v}$, the internal energy ${U_{int}}$ of the bubble content and the potential energy ${E_{pot}}$ of the expanded or collapsed bubble. The potential energy is, at each stage, given by the work done against (or by) hydrostatic pressure, ${W_{stat}}$, plus the work done against (or by) surface tension, ${W_{surf}}$. The dashed arrows represent energy dissipation in each phase via viscous damping, ${W_{visc}}$, vapour condensation, ${E_{cond}}$ and shock wave emission, ${E_{SW}}$. The rebound shock wave emission is driven partly by the internal energy of the collapsed bubble, and partly by the energy stored in the liquid compressed during the collapse phase, $E_{compr}^{coll1}$. Correspondingly, the rebound shock wave energy is composed of two fractions: ${E_{SWB}}$ arising from the bubble rebound, and ${E_{SWL}}$ arising from the re-expansion of the compressed liquid. A complete energy balance can be established only at particular times, when the kinetic energy is zero, i.e. at R_{max 1}, R_{min 1} and R_{max 2}. The energy of the shock wave emitted after optical breakdown is obtained by subtracting the balance established for R = R_{max 1} from ${E_{abs}}$, and the energy of the shock wave emitted during the rebound is evaluated by comparing the balance for R = R_{max 2} with the total energy of the compressed bubble and liquid upon collapse.

The expanding bubble content does work, W_gas, 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

(3.42)\begin{equation}{E_{abs}} = {E_v} + \Delta {U_{int}}(t) + {W_{gas}}(t).\end{equation}

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 ΔU_int into W_gas 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 p_surf arising from the surface tension at the bubble wall. In doing this, it creates kinetic energy E_kin of the accelerated liquid, potential energy E_pot of the expanding bubble, drives the emission of a shock wave with energy E_SW and does the work W_visc by overcoming viscous damping. Altogether, the work done on the liquid involves the components

(3.43)\begin{equation}{W_{gas}} = E_{SW}^{bd} + {E_{kin}} + {W_{visc}} + {E_{pot}},\quad \textrm{with}\;{E_{pot}} = {W_{stat}} + {W_{surf}},\end{equation}

where W_stat and W_surf denote the work done against hydrostatic pressure and surface tension.

At R = R_max, 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 E_abs.

The vaporization energy E_v 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 R_{max 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 R_{max 1} with the amount contained in a bubble with equilibrium radius R_{nc 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 R_{nc 1} and R_{nc 2} = R_nres.

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 E_SWB and E_SWL, 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 E_acoust. 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 R_nres 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 (p_{gas 1}, V ₁) to state 2 (p_{gas 2}, V ₂) is

(3.44)\begin{equation}\Delta {U_{int}} = \frac{{4{\rm \pi}}}{{3(\kappa - 1)}}({p_{gas2}}R_2^3 - {p_{gas1}}R_1^3).\end{equation}

The energy absorbed during plasma formation is deposited into a small volume with equivalent spherical radius R ₀ that is regarded as initial size of the laser-induced cavitation bubble. We assume a pressure

(3.45)\begin{equation}{p_{gas|R = {R_0}}} = {p_\infty } + 2\sigma /{R_0},\end{equation}

for the bubble nucleus at time t ₀ = 0 before the laser pulse, and express the time-dependent gas pressure during the pulse through the time evolution of equilibrium radius R_n as

(3.46)\begin{equation}{p_{gas}}(t) = \left( {{p_\infty } + \frac{{2\sigma }}{{{R_n}(t)}}} \right){\left( {\frac{{R_n^3(t)}}{{{R^3}(t)}}} \right)^\kappa }.\end{equation}

The internal energy of the bubble at time t during bubble expansion is then given by

(3.47)\begin{equation}\Delta U_{int}^{exp}(t) = \frac{{4{\rm \pi}}}{{3(\kappa - 1)}}\left\{ {\left( {{p_\infty } + \frac{{2\sigma }}{{{R_n}(t)}}} \right)\left( {\frac{{R_n^{3\kappa }(t)}}{{{R^{3\kappa }}(t)}}} \right){R^3}(t) - \left( {{p_\infty } + \frac{{2\sigma }}{{{R_0}}}} \right)R_0^3} \right\}.\end{equation}

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 = R_{max 1}. After deleting the second term in (3.47) and considering the van der Waals hard core, it reads

(3.48)\begin{equation}\Delta U_{int}^{coll1}(t) = \frac{{4{\rm \pi}}}{{3(\kappa - 1)}}\left\{ {\left( {{p_\infty } + \frac{{2\sigma }}{{{R_n}(t)}}} \right){{\left( {\frac{{R_n^3(t) - R_{v\,\textrm{d}W}^3}}{{{R^3}(t) - R_{v\,\textrm{d}W}^3}}} \right)}^\kappa } \times ({R^3}(t) - R_{v\,\textrm{d}W}^3)} \right\}.\end{equation}

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

(3.49)\begin{equation}U_{int}^{max1}{|_{p = {p_v}}} = \frac{{4{\kern 1pt} {\rm \pi}}}{{3(\kappa - 1)}}{p_v}R_{max1}^3 = 4{\rm \pi}{p_v}R_{max1}^3.\end{equation}

The internal energy loss during bubble expansion is given by the difference between the energy of an adiabatically expanding bubble at R_max (which is obtained by evaluating (3.47) at the time of maximum bubble expansion for R_n = R_nbd) and the energy corresponding to isothermal conditions at R_max from (3.49)

(3.50)\begin{equation}\Delta U_{int,cond}^{exp} = {U_{int}}{|_{{R_n} = {R_{nbd}}}} - U_{{\mathop{\rm int}} }^{max1}{|_{p = {p_v}}}.\end{equation}

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 R_n = R_{nc 1}) from the energy corresponding to isothermal conditions at R_{max 1}

(3.51)\begin{equation}\Delta U_{int,cond}^{coll1} = U_{int}^{max1}{|_{p = {p_v}}} - {U_{int}}{|_{{R_n} = {R_{nc1}}}}.\end{equation}

The changes of internal energy during rebound and second collapse are

(3.52)\begin{equation}\Delta U_{int,cond}^{reb} = {U_{int}}{|_{{R_n} = {R_{nc1}}}} - U_{int}^{max2}{|_{p = {p_v}}},\end{equation}

and

(3.53)\begin{equation}\Delta U_{int,cond}^{coll2} = \Delta U_{int,cond}^{res} = U_{int}^{max2}{|_{p = {p_v}}} - {U_{int}}{|_{{R_n} = {R_{nc2}}}},\end{equation}

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 R_{max 1}, R_{min 1} and R_{max 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

(3.54)\begin{equation}{\rho _\infty }{\textstyle{4 \over 3}}{\rm \pi} R_0^3 = {\rho _v}{\textstyle{4 \over 3}}{\rm \pi}{(R_v^{bd})^3},\quad \textrm{we}\ \textrm{obtain}\;R_v^{bd} = {R_0}{({\rho _\infty }/{\rho _v})^{1/3}}.\end{equation}

The mass density of vapour at p_v = 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 R_{max 1} and R_{max 2} equals the equilibrium vapour pressure at room temperature, p_v = 2.33 kPa (Lauterborn & Kurz 2010). The corresponding bubble radii for vapour at ambient pressure are then

(3.55)\begin{equation}R_v^{max\,i} = {R_{max\,i}}{({p_v}/{p_\infty })^{1/3}},\quad \textrm{with}\;i = 1\ \textrm{and}\;2.\end{equation}

The loss of latent heat by condensation during bubble expansion is given by

(3.56)\begin{equation}\Delta E_v^{exp} = {E_v} - E_v^{max1},\quad \textrm{with}\;E_v^{max1} = {\left( {\frac{{R_v^{max1}}}{{R_v^{bd}}}} \right)^3}{E_v},\end{equation}

and E_V 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

(3.57a,b)\begin{equation}\Delta E_v^{coll1} \!=\! E_v^{max1} - E_v^{coll1},\quad \textrm{with}\;E_v^{coll1} = {\left( {\frac{{R_v^{coll1}}}{{R_v^{bd}}}} \right)^3}{E_v}\quad \textrm{and}\quad R_v^{coll1} \!=\! {(R_{nc1}^3 - R_{nc2}^3)^{1/3}}.\end{equation}

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

(3.58)\begin{equation}\Delta E_v^{reb} = E_v^{coll1} - E_v^{max2},\quad \textrm{with}\;E_v^{max2} = {\left( {\frac{{R_v^{max2}}}{{R_v^{bd}}}} \right)^3}{E_v},\end{equation}

and

(3.59)\begin{equation}\Delta E_v^{coll2} = E_v^{max2},\end{equation}

respectively. The latent heat remaining at R_{max 2} is dissipated during the transition into linear bubble oscillations. For each stage, the total condensation loss E_cond is the sum of internal energy loss and release of latent heat

(3.60)\begin{equation}{E_{cond}} = \Delta {U_{int,cond}} + \Delta {E_v}.\end{equation}

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

(3.61)\begin{equation}{W_{gas}}(t) = \int {{p_{gas}}(t)\,\textrm{d}V} = \int {4{\rm \pi}\,{R^2}{p_{gas}}(t)\,\textrm{d}R} .\end{equation}

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

(3.62)\begin{equation}{W_{gas}}(t) = \int_0^t {4{\rm \pi} R{{(t^{\prime})}^2}U(t^{\prime}){p_{gas}}(t^{\prime})\,\textrm{d}t^{\prime}} ,\end{equation}

where p_gas is given by (3.46).

The work W_gas done on the liquid during bubble expansion (or by the liquid on the collapsing bubble) is composed of W_stat, W_surf, and W_visc, as expressed by (3.43). The work done at a given time to overcome the hydrostatic pressure is

(3.63)\begin{equation}{W_{stat}}(t) = \int {{p_\infty }\,\textrm{d}V} = \int {4{\rm \pi}{R^2}{p_\infty }\,\textrm{d}R} = \int_0^t {4{\rm \pi}R{{(t^{\prime})}^2}U(t^{\prime}){p_\infty }\,\textrm{d}t^{\prime}} .\end{equation}

The work done to overcome the surface tension is

(3.64)\begin{equation}{W_{surf}}(t) = \int {{p_{surf}}\,\textrm{d}V} = \int {4{\rm \pi}{R^2}{p_{surf}}\,\textrm{d}R} = 8{\rm \pi}\sigma \int_0^t {R(t^{\prime})U(t^{\prime})\,\textrm{d}t^{\prime}} ,\end{equation}

and the potential bubble energy at any given time is therefore

(3.65)\begin{equation}{E_{pot}}(t) = {W_{stat}}(t) + {W_{surf}}(t).\end{equation}

Finally, the work required to overcome viscosity is

(3.66)\begin{equation}{W_{visc}}(t) = \int {{p_{visc}}\,\textrm{d}V} = \int {4{\kern 1pt} {\rm \pi}{R^2}{p_{visc}}\,\textrm{d}R} = 16{\rm \pi}\mu \int_0^t {R(t^{\prime})U{{(t^{\prime})}^2}\,\textrm{d}t^{\prime}} .\end{equation}

The energy balance for the conversion of the absorbed energy E_abs during the expansion process is obtained by numerical integration of the above equations up to the time t_{max 1}

(3.67)\begin{equation}{E_{abs}} = \underline {E_{SW}^{bd}} + \underline {E_{cond}^{exp\,}} + \underline {{W_{visc}}} + E_{pot}^{max1} + U_{int}^{max1}{|_{p = {p_v}}} + E_v^{max1}.\end{equation}

The underscored terms are energy parts that are dissipated during expansion, the other terms remaining at t = t_{max 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 R_{min 1} and tracks their dissipation and conversion during collapse. It is obtained by numerical integration of (3.58)–(3.62) up to t = t_{min 1}

(3.68)\begin{align}E_B^{max1} = E_{pot}^{max1} + U_{int}^{max1}{|_{p = {p_v}}} + E_v^{max1} = \underline {E_{cond}^{coll1}} + \underline {{W_{visc}}} + E_{compr}^{coll1} + U_{int}^{min1} + E_v^{coll1}.\end{align}

The underscored terms represent energy parts that are dissipated during collapse, the other terms remaining at t = t_{min 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 R_{max 2} is determined by numerical integration up to t = t_{max 2}

(3.69)\begin{align}E_{compr}^{total} = E_{compr}^{coll1} + U_{int}^{min1} + E_v^{coll1} = \underline {E_{SW}^{reb}} \!+\! \underline {E_{cond}^{reb}} + \underline {{W_{visc}}} \!+\! E_{pot}^{max2} \!+\! U_{int}^{max2}{|_{p = {p_v}}} + E_v^{max2}.\end{align}

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,

(3.70a,b)\begin{equation}{E_{SWL}} = E_{compr}^{coll1}\quad \textrm{and}\quad {E_{SWB}} = E_{SW}^{reb} - {E_{SWL}}.\end{equation}

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

(3.71)\begin{equation}E_B^{max2} = E_{pot}^{max2} + U_{int}^{max2}{|_{p = {p_v}}} + E_v^{max2} = \underline {{E_{acoust}}} + \underline {E_{cond}^{coll2}} + \underline {{W_{visc}}} + U_{int}^{res}{|_{R = {R_{nc2}}}}.\end{equation}

The energy E_acoust 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 R_{max 3} to determine R_{nc 2}. For nanobubbles, experimental values of R_{max 3} may not be available. In that case, we identify the equilibrium radius during afterbounces and the residual bubble radius with R_{nc 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 ₀ 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 E_th,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 × E_th,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 ₀ for the numerical simulations. Figure 4(d) shows the pulse energy dependence of plasma transmittance T_opt, from which the absorbed energy is determined as ${E_{abs}} = {E_L}(1 - {T_{tra}})$.

Figure 4. Determination of the plasma size produced with 350 fs, 1040 nm pulses of different energy focused at NA = 0.8. (a) Photographs of plasma luminescence taken with the set-up of figure 1(a) and integrated over 70 laser pulses at ISO 3200. (b) Plasma volume determined from photographs as a function of laser pulse energy, assuming rotational symmetry of the plasma around the laser beam axis. (c) Energy dependence of the radius of a sphere having the same volume as the plasma in (b). (d) Plasma transmission.

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.

Figure 5. Initial phase of shock wave emission and bubble expansion after plasma formation produced by 1064 nm, 6 ns laser pulses with pulse energies of (a) 10 mJ and (b) 20 mJ that were focused at NA = 0.25. The laser light was incident from the right. Scale bars represent 100 μm. The self-luminescent plasma appears on all images because photographs were taken with open shutter in a darkened room. Bright-field illumination was done with a 6 ns laser pulse from a collimated laser beam as shown in figure 1(c).

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 2002). 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 1957). 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. 1988; Vogel et al. 1996b). 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. 1996b). 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 1957), 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. 1996a) and streak photography (Noack & Vogel 1998). 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. (1996a). 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 1996; Vogel & Venugopalan 2003; Vogel et al. 2005), the instable liquid starts to expand immediately (Zhigilei et al. 2003). 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.

Figure 6. Confocal probe beam forward scattering signal from a bubble with 35.8 μm maximum radius produced by a 265 fs, 755 nm laser pulse of 155 nJ energy focused at NA = 0.9 The dimensionless stand-off distance from the microscope objective's front lens is γ = 70. (a) Entire signal portraying the transition from nonlinear cavitation bubble oscillations to linear oscillations of the residual gas bubble. (b) Enlarged views of the first four oscillations and the part from 30th to 52nd oscillation. The signal undulations during the first oscillation are interference fringes reflecting the radius–time evolution. The arrows mark the time interval around R_max. Later, each undulation represents one period of the small-amplitude oscillations of the residual bubble. In (c), the oscillation time, T_osc, is plotted as a function of the oscillation number, i. The experimental data are fitted with an asymptotic regression model curve given by ${T_{osc}} = a - b \times {c^i}$, with fitting coefficients a = 0.817, b = −42.267 and c = 0.134. The mean oscillation time from 50th to 102nd oscillation is 813.3±6.7 ns.

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 ₀ = 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⁻³. The internal energy density ${U_{int}}/{V_p} = ({E_{abs}}-{E_v})/{V_p}$ with E_v = 25.49 nJ from (3.41) is 6.14 kJ cm⁻³, and the difference of both values corresponds to the vaporization enthalpy. Since energy deposition is isochoric, the average plasma temperature can be determined from U_int/V_p and the water EOS. Using the IAPWS-95 formulation (Wagner & Pruss 2002), we obtain T_avg = 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 T_{osc 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.

Figure 7. Time evolution of (a) bubble radius, (b) internal pressure and (c) bubble wall velocity corresponding to the signal of figure 6. The inset in (a) shows the laser pulse shape assumed in the calculations and the increase of R_n during the pulse from R ₀ to R_nbd. The reduction of equilibrium bubble pressure at R = R_max shown in (a) goes along with a drop of internal bubble pressure that represents the net amount of vapour condensation during the first bubble oscillation. Peak pressures upon breakdown and collapse and peak velocities are indicated in the figure.

Figure 8. Enlarged views of the time evolution of bubble radius shown in (a) and (b), internal pressure shown in (c) and (d) and bubble wall velocity shown in (e) and (f) after breakdown and around the first bubble collapse for the same parameters as in figure 7. The displayed time interval is 20 ns for the bubble growth and 1 ns for the collapse–rebound phase. The dashed line in the R(t) plot for the collapse phase in (b) represents the van der Waals hard core. Insets in (c) and (e) show the time evolution of P(t) and U(t) during the laser pulse.

Table 1. Characteristic bubble parameters corresponding to the probe beam signal of figure 6 that were used in the simulations of figures 7–10. The equilibrium radii R_n were determined by fitting the predicted R(t) curve to the measured oscillation periods. Specifically, R_nbd was used for fitting T_{osc 1}, R_{nc 1} for fitting T_{osc 2}, and R_{nc 2} for T_{osc 3}. Afterwards, R_n was kept constant for the rest of the calculation. The vapour bubble radius after breakdown corresponds to the vaporized liquid volume with radius R ₀ and is given by (3.52). At R_maxi, it represents the amount of vapour contained in the expanded bubble at the vapour pressure under ambient condition, ${p_v} = 2.33\ \textrm{kPa}$, and is determined using (3.53). At first collapse, it is calculated as ${R_{vc1}} = {(R_{nc1}^3 - R_{nc2}^3)^{1/3}}$, assuming that the bubble content at the second collapse consists only of non-condensable gas. The respective values for the vapour mass m_v were calculated with ${\rho _v} = 0.761\ \textrm{kg}\;{\textrm{m}^{ - 3}}$.

Figure 9. Shock wave emission after breakdown for the parameters of figure 7, with velocity distributions in the liquid, u(r), at different times in (a), and the corresponding pressure distributions, p(r), presented in (b). The circles indicate the respective velocity and pressure values at the bubble wall and its position. The inset in (a) shows an enlarged view of the shock wave propagation when it has detached from the outward going radial flow in the bubble's vicinity. The dash-dotted line in (b) represents a decay curve of the shock wave's peak pressure, p_peak(r) that was derived from 144 p(r) profiles. The slopes of the p_peak(r) curve are indicated for various propagation distances. The pressure decay is faster than for acoustic waves for which the attenuation would be proportional to r ⁻¹.

Figure 10. (a) Evolution of the velocity distribution in the liquid during the late stage of bubble collapse and during the bubble's rebound for the parameters of figure 7. The time evolution of the u(r) curves is shown with the circles representing the respective pressures at the bubble wall and its position. The times given for the individual u(r) curves refer to the instant at which the bubble reaches its minimum radius, which is set as t = 0. On a time scale starting with bubble generation, it corresponds to t_coll = 6488.0 ns. Upon rebound, the flow around the expanding bubble collides with the still incoming flow from outer regions, and a shock front develops within about 50 ps and 750 nm propagation distance that continues to exist even in the far field. (b) Evolution of the pressure distribution in the liquid. The times given for the individual p(r) curves refer to the instant at which the bubble reaches its minimum radius, which is set as t = 0. The curve for bubble wall position during the collapse phase was determined from 28 shock wave profiles, and the respective curve for the rebound phase as well as the p_peak(r) curve were derived from 53 p(r) profiles. After the shock front has formed, the amplitude of the outgoing pressure wave drops initially very rapidly and later more slowly. However, even in the far field, the shock front continues to exist and the pressure decay is faster than for acoustic waves.

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. 1996a). Water dissociation relies on free-electron-mediated pathways as thermal dissociation sets in only at 3000 K (Mattsson & Desjarlais 2006). 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 R_min = 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 2011; Koch et al. 2016). 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⁻¹ (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 = R_max 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 ₀. 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⁻¹ to a peak rebound value of 369 m s⁻¹ (figure 8f) which goes along with an acceleration of 5.4 × 10¹³ m s⁻². 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 (1957) 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⁻¹ 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 R_{nc 1} is T_coll = 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 p_coll = 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 ⁻¹. 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 1954; Rogers 1977).

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. 1999b). 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 1957). 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 R_v after plasma formation is larger than R_nbd. The entity R_nbd 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 R_v and R_nbd 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 R_{vc 1} is not much smaller than R_{nc 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 R_nc 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 2000). However, during second collapse almost all vapour will condense because the rebounded bubble at R_{max 2} contains much less vapour than the larger bubble at R_{max 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 R_{nc 2} = 2.44 μm, while R_{nc 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, R_res, 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 R_nres = 2.60 μm for the equilibrium radius of the residual gas bubble under isothermal conditions (κ = 1). Note that R_nres = 2.60 μm, is slightly larger than R_{nc 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, T_W, 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.

Figure 11. Simulated R(t) curves for the signal of figure 6 covering the transition from nonlinear to linear oscillations and late bubble oscillations. (a) Calculations with values for surface tension σ and viscosity μ at room temperature, T_W = 20 °C. The inset shows an expanded view around the equilibrium gas bubble radius R_res for times between 25 and 35 μs. Simulations were performed for adiabatic conditions during the entire bubble lifetime, with κ = 4/3, and for initially adiabatic conditions followed by isothermal conditions after the maximum of the third oscillation, with κ = 1. (b) Shows R(t) curves for T_W = 60 °C and 110 °C, with room temperature as reference. Simulations were performed for initially adiabatic conditions followed by isothermal conditions after the maximum of the third oscillation. The inset shows an expanded view around R_res for times between 40 and 70 μs. (c) Shows the R(t) curve for times up to the end of the experimentally observed bubble oscillations at 90 μs, with further expanded radius scale.

The simulation results for room temperature (T_W = 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 T_W = 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 R_res = 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 E_abs = 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 T_avg = 68.7 °C is consistent with a value T_W = 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 T_{osc, 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 (R_{max 1}, R_{min 1} and R_{max 2}), it is possible to establish a complete energy balance because the fluid movement has stopped, E_kin = 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.

Table 2. Energy balance for the signal of figure 6. Terms in brackets denote subfractions of the energy parts listed above them. Laser and bubble parameters are the same as in figures 7–10. The table lists all energy values contained in the flow diagram of figure 3. Starting point is the absorbed laser energy, which at the end of the laser pulse partitions into ${E_{abs}} = {E_v} + \Delta {U_{int}}$, with E_v = 25.49 nJ (40.8 %) and ΔU_int = 37.0 nJ (59.2 %). The first column shows how E_abs set as 100 % splits into various parts, whereby the fractions dissipated during expansion are listed in the upper part of the column, and the parts remaining at R = R_{max 1} are given in the lower part. The energy of the expanded bubble, $E_B^{max1}$, is the starting point for partitioning in the collapse phase, whereas the energy $E_{compr}^{total}$ of the compressed bubble content and compressed liquid at R = R_{min 1} is the reference point for the rebound phase, and $E_B^{max2}$ for the afterbounces.

Figure 12. Time evolution of the energy partitioning of internal bubble energy U_int for the same parameters as in table 2. Panel (a) shows the increase of work W_gas done on the liquid by the expanding gas together with the corresponding decrease of internal energy, followed by a reversal upon collapse, where the inrushing liquid does work on the bubble content. The inset shows a magnified view of the collapse and of the evolution of bubble energy up to the fifth oscillation. The shock wave energy is included in W_gas but can be explicitly evaluated only when E_kin = 0 at R_{max 1} and R_{max 2} (see table 2). Panel (b) presents on an expanded scale the evolution of the energy fractions needed to overcome viscous damping, W_visc, and surface tension, W_visc, as well as the internal energy lost via condensation, E_cond. This change of E_cond reflects the changes of the bubble's equilibrium radius during the bubble oscillations from R_nbd through R_{nc 1} to R_{nc 2}. In reality, these changes reflecting condensation occur continuously, and the jumps at R_{max 1} and R_{max 2} in (b) are due to the simplified way in which they are considered in the present model.

The incident energy of the laser pulse producing the bubble with R_max = 35.8 μm was 155 nJ, and the absorbed energy determined from the measured plasma transmission was E_abs = 81 nJ (§ 4.1). The total amount of deposited energy derived from the plasma radius R ₀ 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 E_abs) appears as potential energy $E_{pot}^{max1}$ of the expanded cavitation bubble at R_{max 1}, and 38.9 % of $\Delta {U_{int}}$ (of E_abs) 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 R_min. 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. 2001; Byun & Kwak 2004; Dagallier et al. 2017). 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 ₀ is obtained by identifying its volume with the breakdown volume, and model predictions for T_osci are fitted to experimental values by varying the equilibrium bubble radii R_nbd, R_{nc 1} and R_{nc 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 R_n. 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. (1996a) 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⁻³, 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.

Figure 13. Comparison of simulated U(t) curves after breakdown with the experimental data from Vogel et al. (1996b) for laser-induced bubbles generated at 1064 nm wavelength. Pulse durations and energies were 30 ps and 50 μJ in (a), 30 ps and 1 mJ in (b), 6 ns and 1 mJ in (c) and 6 ns and 10 mJ in (d). Simulations were performed with and without consideration of the contribution of particle velocity behind the shock front to the bubble wall velocity. Simulation parameters providing an optimum fit to experimentally determined R_{max 1} values are R ₀ = 8.5 μm, R_nbd = 86.1 μm, R_nbd = 87.2 μm in (a), R ₀ = 26 μm, R_nbd = 294 μm, R_nbd = 298.3 μm in (b), R ₀ = 19 μm, R_nbd = 291 μm, R_nbd = 297 μm in (c) and R ₀ = 37 μm, R_nbd = 660 μm and R_nbd = 671 μm in (d).

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⁻¹ using (3.11) that is based on Hugoniot data valid up to 25 GPa (Rice & Walsh 1957). The Tait equation (3.4) used in the Gilmore model leads to the relationship

(5.1)\begin{equation}{p_s} = ({p_\infty } + B)\left( {\frac{{2nu_s^2}}{{(n + 1)c_0^2}} - \frac{{n - 1}}{{n + 1}}} \right) - B,\end{equation}

between velocity and pressure at the shock front (Müller 1987; Vogel et al. 1996a). For u_s = 4500 m s⁻¹, (3.11) yields a 1.63 times larger pressure than (5.1). If we correct the Gilmore model predictions by that factor, we obtain p_s = 7.73 GPa, in very good agreement with the experimental value of 7.15 GPa

Figure 14. (a) Simulated P(t) curves and R(t) curves (inset) for the initial phase of laser-induced bubble expansion for the same parameters as in figure 13(d). Solid lines show the results with jump start of the bubble wall velocity, and dashed lines show results without its consideration. (b) Pressure distributions in the liquid at different time instants showing the formation of a shock front and the initial phase of shock wave emission. The circles indicate the pressure P at the bubble wall and its position for the respective p(r) curves, and the dotted lines show the P(R) trajectories. The rapid start of the bubble motion with jump start of the bubble wall velocity results in a lower maximum bubble pressure (4749 MPa) than without jump start (8803 MPa). However, the shock front has formed after about 8 ns in both cases, as seen in (b).

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 1980; Storey & Szeri 2000; Akhatov et al. 2001; Zein et al. 2013; Zhong et al. 2020; Trummler, Schmidt & Adams 2021). The specific nature of the non-condensable gas will differ, depending on the type of cavitation. Akhatov et al. (2001) 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. 1983; Sato et al. 2013; Barmina et al. 2016, 2017).

Water dissociation can proceed through thermally driven hydrolysis (Mattsson & Desjarlais 2006, 2007) and by free-electron-mediated bond breaking. Liang, Zhang & Vogel (2019) 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. 1996; Fedor et al. 2006; Ram, Prabhudesai & Krishnakumar 2009). 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. 2015). The hydrated electrons may also contribute to dissociation and subsequent gas formation (Draganic & Draganic 1971; Nikogosyan et al. 1983; Elles et al. 2007) because their energy (≥6.5 eV) is larger than the O–H bonding energy of § 5.2 eV (Toegel, Hilgenfeldt & Lohse 2002; Maksyutenko, Rizzo & Boyarkin 2006). 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 1983; Mattsson & Desjarlais 2006). For T > 4500 K, most of the chemical bonds are broken, and only radicals of H and O are present (Lédé et al. 1983; Jung, Jang & You 2013). Correspondingly, the amount of non-condensable gas in laser-produced cavitation bubble was found to increase with growing plasma temperatures (Sato et al. 2013).

In the case of laser-induced bubble generation that is portrayed in detail in the present paper, the plasma temperature is relatively low (T_avg = 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. 2014). 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 2000). 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 R_max, 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 2000). 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 1980; Yasui 1995; Storey & Szeri 2000; Toegel et al. 2000; Akhatov et al. 2001; Lauer et al. 2012; Zein et al. 2013; Zhong et al. 2020). Unfortunately, the values of evaporation and condensation coefficients depend on pressure and temperature, and their values are still uncertain (Eames, Marr & Sabir 1997; Marek & Straub 2001). We escaped the dilemma by using R_n 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. (2001).

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⁻³ to approximately 40 kJ cm⁻³, depending on laser pulse duration, wavelength, focusing angle and pulse energy (Vogel et al. 1996b). 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 R_{max 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 × 10⁵. 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. 1999, 2001; Ohl et al. 1999; Brenner et al. 2002; Mattsson & Desjarlais 2006).

Upon collapse, most energy is carried away by shock wave emission and little energy remains for the rebounding bubble, which results in a large R_{max 1}/R_{max 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 R_{max 1}/R_{max 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. 2000). The parameter combination explored in the present study (small bubble and relatively low plasma energy density) results in a particularly strong collapse.

Figure 15. Ratio R_{max 1}/R_{max 2} under ambient conditions as a function of maximum bubble radius. ■ Present study; ▰ Vogel & Lauterborn (1988); □ Akhatov et al. (2001); ● Sinibaldi et al. (2019). The ratio is a measure of energy dissipation during the first collapse. Since most energy is carried away by acoustic radiation, it is also indicative of the amplitude of the collapse pressure. An arrow marks the data point corresponding to the signal of figure 6. The collapse of the highly spherical laser-induced bubbles investigated in this paper is more vigorous than that of larger, millimetre-sized bubbles, where usually smaller focusing angles were used for plasma generation and a combination of elongated plasma shape and buoyancy led to deviations from spherical shape. The largest R_{max 1}/R_{max 2} value for millimetre-sized bubbles was observed by Sinibaldi et al. (2019) for tight focusing at NA = 0.6. The R_{max 1}/R_{max 2} ratio increases for R_{max 1} < 10 μm, where surface tension and viscosity become important.

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. 2000) and with a driving pressure p_a ≈ 0.145 MPa (Matula 1999; Toegel & Lohse 2003). Simulations by Matula (1999) predicted a ratio R_max/R_nc = 9.94 for a SBSL bubble driven at f ≈ 25 kHz and p_a ≈ 0.142 MPa. Since a similar ratio (R_{max 1}/R_{nc 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 (2000) found a shock wave velocity of u_s = 4000 m s⁻¹ at f ≈ 20 kHz and p_a ≈ 0.139 MPa, and Weninger et al. (2000) reported u_s = 5930 m s⁻¹ (Mach 4) for f ≈ 16.5 kHz and p_a ≈ 0.145 MPa. An evaluation of these u_s data using the Hugoniot data of Rice & Walsh (1957) yields p_s = 5.3 GPa and p_s = 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 1996; Lin & Szeri 2001; Brenner et al. 2002). 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 10⁶ K (Wu & Roberts 1993). However, Storey & Szeri (2000) 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. 2001) 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 (1957) is 5400 m s⁻¹. This value is much larger than the peak bubble wall velocity of 1793 m s⁻¹, 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, E_SW during the entire bubble life is 54.2 % of E_abs, and E_cond 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 1957) 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 (1964) and Akulichev et al. (1968) as well as the simulations by Fujikawa & Akamatsu (1980) were based on assumptions about the bubble's gas content rather than on measured R(t) curves or oscillation times. Akhatov et al. (2001) 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 1964; Fujikawa & Akamatsu 1980; Koch et al. 2016). We observe a rapid formation of shock waves as previously reported by Akulichev et al. (1968) and Ebeling (1978), 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. 1996a).

Experimentally observed pressure decays are usually steeper than predictions by the Gilmore model, with steepest slopes ≥2.0 (Vogel et al. 1996a; Lai et al. 2021). 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. (1996a) 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 1988). However, the shock wave propagation in the far field exhibits radial symmetry even when the breakdown region is elongated (Tagawa et al. 2016). 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. 1999b; Tinguely et al. 2012). 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 E_SWL/E_SWB = 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. 2011). 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. (2011) was only 1.5 GPa according to the Gilmore model, and 2.6 GPa according to their full model. The ratio E_SWL/E_SWB 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, E_SWL/E_SWB 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. 1999b). Furthermore, an ever-larger part of the mechanical energy appears as shock wave energy (Lai et al. 2021). For a 6 ns, 10 mJ pulse with ε = 40 kJ cm⁻³, the shock wave energy was found to be more than two times larger than the bubble energy (Vogel et al. 1999b), different from the case of figure 6 in the present paper, where for ε = 8.7 kJ cm⁻³ 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 R_{max 1}/R_{max 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$.