1. Introduction
Encapsulated fluid structures (Sattari, Hanafizadeh & Hoorfar Reference Sattari, Hanafizadeh and Hoorfar2020) are complex bodies made of multiple fluids, where one is enclosed in another. Such structures are commonly classified into three families, gas in liquid, liquid in liquid and liquid in gas, where the outer fluid serves as a ‘shell’. They may also be divided into two categories based on the size ratio between the shell and the encapsulated body: thin-shelled (
$1-\varepsilon <0.002$) and thick-shelled (
$1-\varepsilon >0.2$) structures, where 
$\varepsilon$ is the drop-to-bubble radius ratio. Examples of gas-in-liquid thin- and thick-shelled structures are soap bubbles and drop-encapsulated bubbles (Bartolo, Josserand & Bonn Reference Bartolo, Josserand and Bonn2006), respectively. Liquid-in-liquid structures having a thin and thick shell are named liquid-coated drops (Polenz, Weitz & Baret Reference Polenz, Weitz and Baret2015) or emulsions (Kim & Kim Reference Kim and Kim2014; Lim et al. Reference Lim, Koh, Ren, Chin, Shi and Yan2017), respectively. Gas encapsulating a liquid is generally called an antibubble (Dorbolo, Caps & Vandewalle Reference Dorbolo, Caps and Vandewalle2003; Scheid et al. Reference Scheid, Dorbolo, Arriaga and Rio2012; Vitry et al. Reference Vitry, Dorbolo, Vermant and Scheid2019) when the gas layer is very thin and, more rarely, a gas-encapsulated drop (GED) (Shen et al. Reference Shen, Hu, Chen, Xie and Fu2018; Biasiori-Poulanges et al. Reference Biasiori-Poulanges, Bokman, Baumann and Supponen2022) for larger gas shells. Encapsulated structures have been of particular interest in the fields of fluid transport (Gref et al. Reference Gref, Minamitake, Peracchia, Trubetskoy, Torchilin and Langer1994; Zia et al. Reference Zia, Nazir, Poortinga and van Nostrum2022), chemical microreaction (Song, Chen & Ismagilov Reference Song, Chen and Ismagilov2006), mixing processes (Wang et al. Reference Wang, Jiao, Huang, Yang and Nguyen2009), biomedical applications (Yang et al. Reference Yang, Guo, Kiraly, Mao, Lu, Leong and Huang2012) and targeted drug delivery (Duncanson et al. Reference Duncanson, Abbaspourrad, Shum, Kim, Adams and Weitz2012; Boissenot et al. Reference Boissenot, Bordat, Fattal and Tsapis2016). They are also considered as adequate microcarriers in applications involving the encapsulation of cells (Choi et al. Reference Choi, Wang, Lee, Kim, Zhang, Mao, Mooney and Weitz2016b; Lee et al. Reference Lee, Praveenkumar, Oh, Lee and Kim2016), photonic ink capsules (Kim et al. Reference Kim, Lee, Yang and Yi2011; Shirk et al. Reference Shirk, Steiner, Kim, Marquez and Martinez2013) or the encapsulation of hydrophobic cargo in polymeric microcapsules (Choi et al. Reference Choi, Lee, Abbaspourrad, Kim, Fan, Caggioni, Wesner, Zhu and Weitz2016a).
In the field of targeted drug delivery, microbubbles have often been promoted as one of the most promising candidates for drug transport and release, leveraging ultrasound to induce the payload discharge (Postema & Gilja Reference Postema and Gilja2007; Tinkov et al. Reference Tinkov, Bekeredjian, Winter and Coester2009; Shakya et al. Reference Shakya, Cattaneo, Guerriero, Prasanna, Fiorini and Supponen2024). They have first been clinically approved as ultrasound contrast agents (Kiessling et al. Reference Kiessling, Fokong, Koczera, Lederle and Lammers2012), where a coating typically made of phospholipids or proteins is employed to increase their stability and lifetime. Two main techniques have been proposed to load these vehicles with drugs, either by attaching loaded liposomes (McLaughlan et al. Reference McLaughlan, Harput, Abou-Saleh, Peyman, Evans and Freear2017) or nanoparticles (Mørch et al. Reference Mørch2015) to the exterior of the microbubble shell, or by adding an inner loaded layer of phospholipids or oil to the bubbles. However, these methods suffer from a low-drug-loading capacity, unstable drug discharge and are often limited to hydrophobic drugs (Fokong et al. Reference Fokong2012). A possible answer to these shortcomings are antibubbles, which are currently being promoted as an alternative to microbubbles in the context of targeted drug delivery (Kotopoulis et al. Reference Kotopoulis2022). These differ from the ephemeral antibubbles first described as inverse bubbles (Hughes & Hughes Reference Hughes and Hughes1932; Stong Reference Stong1974) as well as the gas-encapsulated structures (Shen et al. Reference Shen, Hu, Chen, Xie and Fu2018) in which the relative motion of the smaller liquid drop within the gas layer increases its stability against full coalescence up to tens of seconds. The micrometric antibubbles potentially applicable in targeted drug delivery (Poortinga Reference Poortinga2013; Silpe et al. Reference Silpe, Nunes, Poortinga and Stone2013) are stabilised against dissolution by adding hydrophobised silica nanoparticles to their surface (Poortinga Reference Poortinga2011), also called Pickering stabilisation, which increases the lifetime of these liquid structures up to tens of hours. Although they are thick-shelled bodies, Poortinga (Reference Poortinga2011) refer to them as antibubbles and no distinction is made with GEDs. In addition, they often enclose multiple drop cores of variable sizes. The main advantages of encapsulated structures in this context are the increased loading potential, isolation of the core material from the surrounding medium and the possibility of mixing incompatible compounds (Mishra Reference Mishra2015).
Controlled and/or targeted release in drug delivery trials have been performed using ultrasound as a driving mechanism. In biomedical applications, acoustic pressures in the tens or hundreds of kilopascals and a frequency in the megahertz range are typically used to excite microbubbles into motion that transiently increases the permeability of cells and tissue to drugs without negative bioeffects (McDannold, Vykhodtseva & Hynynen Reference McDannold, Vykhodtseva and Hynynen2008; Carpentier et al. Reference Carpentier2016; Moreno-Gomez et al. Reference Moreno-Gomez, Athanassiadis, Poortinga and Fischer2023). Sonopermeabilisation in the context of transdermal or blood–brain barrier opening has also been promoted using shock waves (Ohl et al. Reference Ohl, Arora, Ikink, De Jong, Versluis, Delius and Lohse2006). In fact, shock waves are already in clinical use such as in lithotripsy for treating kidney stones and gall stones or in low-amplitude shock wave therapy for treating tendinopathies (Rompe et al. Reference Rompe, Hopf, Küllmer, Heine, Bürger and Nafe1996). Shock waves have the advantage of triggering an instantaneous and strong response from a bubble most often in the form of a high-speed liquid jet travelling in the direction of propagation of the shock wave (Philipp et al. Reference Philipp, Delius, Scheffczyk, Vogel and Lauterborn1993; Kodama & Takayama Reference Kodama and Takayama1998; Ohl & Ikink Reference Ohl and Ikink2003; Wolfrum et al. Reference Wolfrum, Kurz, Mettin and Lauterborn2003; Bokman et al. Reference Bokman, Biasiori-Poulanges, Meyer and Supponen2023b), which could be beneficial in the context of targeted drug delivery.
Recent investigations of the interaction of micrometric antibubbles with ultrasound have demonstrated a slow release of the drop using pressure amplitudes and frequencies in the range 1–100 kPa and 0.09–1 MHz, respectively (Moreno-Gomez et al. Reference Moreno-Gomez, Athanassiadis, Poortinga and Fischer2023). The use of shock waves has been suggested in a recent numerical study (Biasiori-Poulanges et al. Reference Biasiori-Poulanges, Bokman, Baumann and Supponen2022), potentially offering a faster and more-controlled release of the drop compared to ultrasound.
The present study experimentally and numerically examines the impulsive release of GEDs with a focus on their use in the field of fluidic transport and mixing where a controlled release of the drop is desirable. The interaction of thin shock waves with large GEDs is studied using synchrotron X-ray phase contrast imaging uniquely allowing visual access to the drop dynamics within the bubble during the release process. Numerical simulations complement the experimental observations by providing additional information on the flow fields and drop dispersion upon release. This work aims at identifying and characterising the different shock-wave-driven drop breakup regimes and the subsequent drop release dynamics.
2. Method
2.1. Experimental method
The interaction between GEDs and shock waves in water is investigated using the experimental set-up illustrated in figure 1(a). The multilayer encapsulation process of liquid drops in air is carried out following the procedure described by Shen et al. (Reference Shen, Hu, Chen, Xie and Fu2018), which consists of generating a steady Taylor flow using a custom-made microfluidic device composed of a PEEK low-pressure Tee with a 1-mm through hole (P-714, IDEX), to which a glass capillary of inner diameter 
$d_{c,i} = 0.7$ mm and length 
$l_{c} = 100$ mm is connected to enable the discharge of the GEDs into the water tank. The Tee is positioned vertically and air is injected at the lower inlet while water is provided from the side as depicted in figure 1(a). The air flow is provided by an air compressor (Fatmax DST 101/8, 8 bar, Stanley) and controlled by a precision pressure regulator (RP1000-8G-02, CKD). The water flow is controlled by a syringe pump (Pump 33 DDS Dual Drive System, Harvard apparatus). The water flow rate is set at a constant 
$3.5\ {\rm ml}\ {\rm min}^{-1}$ and the air pressure is steadily increased until the pressure regulator reads 
$10$ kPa, resulting in a consistent encapsulation process. A sequence of radiographs showing the formation of a GED is displayed at the bottom of figure 1(a). Once a gas slug exits the microfluidic device, it forms a bubble which is followed by a shorter liquid slug, which holds enough momentum to deform and partially penetrate the bubble. As the liquid column stretches and decelerates inside the bubble, it starts necking and eventually pinches off a drop. This process is repeatable, with the bubbles displaying a mean equivalent initial radius and standard deviation of 
$r_{b,0} = 1.13 \pm 0.06$ mm and the drops of 
$r_{d,0} = 0.39 \pm 0.09$ mm (number of specimen, 
$N=71$). The equivalent bubble radius is computed from the projected area of deformed bubbles. The bubble and drop sizes are fixed and constrained by the geometry of the microfluidic device, and the drop location within the bubble varies randomly through the random timing of the start of the experiment with respect to the drop motion within the bubble. As GEDs are generated at the bottom of the water-filled tank, they rise due to buoyancy and remain stable over a few centimetres. The drop's relative motion with respect to the bubble prevents its total coalescence with the surrounding water.
Figure 1. (a) Side view of the experimental set-up. An image sequence showing the gas encapsulation of a drop, where the snapshots are taken at a frequency of 500 fps, and the temporal shock wave pressure profile recorded by the hydrophone are displayed. The part of the pressure signal superior to the ambient pressure is highlighted in dark. (b) Top view of the X-ray beamline and dual imaging system of the ID19 at the European Synchrotron Radiation Facility.
Shock waves are induced in water by optical breakdown. A 6-mm-diameter pulsed laser beam, originating from a Nd:YAG laser (Q-smart, 532 nm, 220 mJ, 5 ns, Lumibird), is expanded through a 
$10\times$ beam expander and focused into a single point by a 90
$^\circ$ parabolic mirror (Aluminum Off-Axis Mirror, Edmund Optics) ionising the liquid. The resulting plasma emits a shock wave as it expands and develops a vapour bubble. The origin of the shock wave is located 5–7 mm away from the stream of GEDs. The pressure waveform is measured by a 75-
$\mathrm {\mu }$m needle hydrophone (NH0075, Precision Acoustics) located 40.5 mm away from the shock wave origin, and recorded using an oscilloscope (WaveRunner 9000, Teledyne LeCroy). The shock wave profile is a sharp, step-like increase in pressure, followed by an exponentially decaying tail of very short duration (see inset of figure 1a), much weaker than the ones of infinite duration reported in the work of Biasiori-Poulanges et al. (Reference Biasiori-Poulanges, Bokman, Baumann and Supponen2022). The duration during which the shock wave exhibits a pressure superior to the ambient pressure is 
$t_{I}\approx 1\ \mathrm {\mu }$s (see highlighted portion in figure 1a). Such shock waves are referred to as impulsive shock waves (Tomita & Shima Reference Tomita and Shima1986; Bokman et al. Reference Bokman, Biasiori-Poulanges, Meyer and Supponen2023b) because of the few microseconds they act on the GED in contrast to the hundreds of microseconds it takes for the ensuing bubble dynamics to unfold. The peak pressure and pressure impulse are tuned by adjusting the laser energy within ranges of 
$p_{max} = 31.38$–93.23 MPa and 
$j =2.63$–5.44 Pa s, respectively. The pressure impulse is obtained by integrating the pressure waveform of the shock with respect to time over the duration 
$t_{I}$. Further details can be found in Bokman et al. (Reference Bokman, Biasiori-Poulanges, Meyer and Supponen2023b). The maximum Mach number encountered in the experiments is 
$M = u_{shock}/c_{0} = 1.06$, estimated using the Rankine–Hugoniot jump relations (Vogel, Busch & Parlitz Reference Vogel, Busch and Parlitz1996), where 
$u_{shock}$ is the propagation speed of the shock wave and 
$c_{0}= 1481\ {\rm m}\ {\rm s}^{-1}$, the speed of sound in water at 20
$\,^\circ$C. The ambient pressure is 
$p_{0} = 101\ 325$ Pa.
In-situ high-speed synchrotron-based X-ray phase contrast imaging is used to obtain optical access to all phase discontinuities along the X-ray beam path and thereby to the drop inside the bubble. Radiographs are recorded at the European Synchrotron Radiation Facility (ESRF) 150 m-long ID19 beamline. The polychromatic hard X-ray beam (mean energy 30 keV) is generated using two axially aligned long-period undulators and conditioned with a series of filters and in-vacuum slits along the vacuum flight tube, moderating the heat load delivered to the water tank located in the experimental hutch. Figure 1(b) shows the X-ray beam travelling through the 
$100\times 100\times 350\ {\rm mm}^{3}$ custom-made water container specifically designed to reduce the X-ray absorption through the use of telescopic windows (SM1L10, Thorlabs). The telescopic windows can slide along the axis of the X-ray beam and the inter-window distance, 
$d_{w}$, is set to 20 mm, above which the water absorption prevents detection of the transmitted signal. The transmitted beam is converted into visible light utilising the 500-
$\mathrm {\mu }$m-thick LYSO:Ce scintillator. Two pellicle beamsplitters are used to reorient the visible light onto two synchronised ultra-high-speed Shimadzu Hyper Vision HPV-X2 cameras (Cameras 1 and 2). The cameras are equipped with 1
$\times$ magnification (i.e. 
$32\ \mathrm {\mu }{\rm m}\ \mathrm {px}^{-1}$) and are set in series to double the acquisition time at a frequency of 0.568 Mfps and an exposure time of 1460 ns (Escauriza et al. Reference Escauriza2020). The detector arrangement is located 7.5 m downstream of the water container, ensuring sufficient X-ray phase contrast contribution through free-space propagation. An increased edge contrast due to the (partial) spatial coherent illumination is therefore achieved while preserving the bubble and drop shape on the images (Wilkins et al. Reference Wilkins, Gureyev, Gao, Pogany and Stevenson1996). Further information on the X-ray imaging system is provided by Bokman et al. (Reference Bokman, Biasiori-Poulanges, Lukić, Bourquard, Meyer, Rack and Supponen2023a). Note that the relatively small dimensions of the water container and inter-window gap result in the GEDs being affected by reflections of shock waves and expansion waves off boundaries. The first reflections from the windows and container for a pressure discontinuity propagating close to 
$M=1$ are estimated to influence a GED located at the set-up's centre slightly before 14 
$\mathrm {\mu }$s and 68 
$\mathrm {\mu }$s, respectively. However, considering the dissipation of the shock wave energy as it travels in water, the influence of reflections is deemed secondary. The influence of gravity is neglected considering the difference in timescales between the drop motion within the bubble (ms) and drop release mechanism (
$\mathrm {\mu }$s).
2.2. Numerical method
2.2.1. Governing equations
The diffuse interface method applied herein uses a thermodynamically well-posed, pressure- and temperature-disequilibrium, multicomponent flow model conserving mass, momentum and total energy (Schmidmayer et al. Reference Schmidmayer, Petitpas, Daniel, Favrie and Gavrilyuk2017). The model and its associated numerical method are implemented in the open-source code ECOGEN (Schmidmayer et al. Reference Schmidmayer, Petitpas, Le Martelot and Daniel2020), which have been validated for several compressible multiphase flow configurations, including bubbles (Pishchalnikov et al. Reference Pishchalnikov, Behnke-Parks, Schmidmayer, Maeda, Colonius, Kenny and Laser2019; Trummler et al. Reference Trummler, Bryngelson, Schmidmayer, Schmidt, Colonius and Adams2020) and drops (Dorschner et al. Reference Dorschner, Biasiori-Poulanges, Schmidmayer, El-Rabii and Colonius2020). The different fluid components are assumed to be immiscible. Although capillary and viscous effects are three to four orders of magnitude lower than inertial effects in the shock-induced dynamics of the bubble and drop, surface tension is activated as it helps stabilising spurious shock-driven interfacial instabilities on the bubble. The model reads
\begin{equation} \left.\begin{array}{c@{}}
\dfrac{\partial \alpha_{i}}{\partial t} + \boldsymbol{u}
\boldsymbol{\cdot} \boldsymbol{\nabla} \alpha_{i} =
\delta p_{i}, \\ \dfrac{\partial \alpha_{i}
\rho_{i}}{\partial t} + \boldsymbol{\nabla}
\boldsymbol{\cdot} ( \alpha_{i} \rho_{i}
\boldsymbol{u}) = 0, \\ \dfrac{\partial \rho
\boldsymbol{u}}{\partial t} + \boldsymbol{\nabla}
\boldsymbol{\cdot} ( \rho \boldsymbol{u} \times
\boldsymbol{u} + p \boldsymbol{I} + \boldsymbol{\varOmega})
= \boldsymbol{0}, \\ \dfrac{\partial \alpha_{i} \rho_{i}
e_{i}}{\partial t} + \boldsymbol{\nabla}
\boldsymbol{\cdot} ( \alpha_{i} \rho_{i} e_{i}
\boldsymbol{u} ) + \alpha_{i} p_{i} \boldsymbol{\nabla}
\boldsymbol{\cdot} \boldsymbol{u} ={-}p_{I} \delta
p_{i} , \\ \dfrac{\partial f}{\partial t} + \boldsymbol{u}
\boldsymbol{\cdot} \boldsymbol{\nabla} f = 0 ,
\end{array}\right\} \end{equation}
where 
$\alpha _{i}$ is the volume fraction, 
$\boldsymbol {u}$ is the flow velocity vector, 
$p_{i}$ is the pressure, 
$\rho _{i}$ the density, 
$e_{i}$ is the internal energy and 
$\delta p_{i}$ is the pressure-relaxation term of the 
$i$th phase. For an 
$N$-phase mixture, the density, pressure and mixture total energy are
\begin{equation} \left.\begin{array}{c@{}}
\displaystyle \rho = \sum_{i=1}^{N} \alpha_{i}\rho_{i}, \\
\displaystyle p =\sum_{i=1}^{N} \alpha_{i}p_{i}, \\
\displaystyle E = \sum_{i=1}^{N} Y_{i}e_{i}(\rho_{i},p_{i})
+ \dfrac{1}{2}\|\boldsymbol{u}\|^{2} +
\dfrac{\varepsilon_\gamma}{\rho}, \end{array}\right\}
\end{equation}
respectively. Here, the first term of the mixture total energy is the mixture-specific internal energy, 
$Y_{i} = \alpha _{i}\rho _{i}/\rho$ is the mass fraction of phase 
$i$ and its internal energy, 
$e_{i}(\rho _{i},p_{i})$, is defined using an equation of state. The last term is the capillary energy 
$\varepsilon _\gamma = \gamma \| \boldsymbol {\nabla } f \|$, where 
$\gamma$ is the surface-tension coefficient, and 
$f$ is a colour function.
The capillary tensor reads
\begin{equation} \boldsymbol{\varOmega} ={-} \gamma \left( \| \boldsymbol{\nabla} f \| \boldsymbol{I} - \frac{\boldsymbol{\nabla} f \times \boldsymbol{\nabla} f}{\| {\boldsymbol{\nabla} f} \|} \right) . \end{equation}GEDs interacting with laser-induced shock waves are investigated using three phases, modelled as a mixture of air, vapour and water for which the air and vapour are described using the ideal-gas equation of state
\begin{equation} p_{i} = \rho_{i}(\kappa_{i}-1)(e_{i}-e_{i,ref}), \end{equation}
where 
$\kappa _{i}$ and 
$e_{i,ref}$ are model parameters corresponding to 
$\kappa _{a} = 1.4$ and 
$e_{a,ref} = 0\ \textrm {J}\ \textrm {kg}^{-1}$ for air and 
$\kappa _{v} = 1.43$ and 
$e_{v,ref} = 2.03\ \textrm {MJ}\ \textrm {kg}^{-1}$ for vapour (Le Métayer, Massoni & Saurel Reference Le Métayer, Massoni and Saurel2004). The reference densities are set to 
$\rho _{a} = 1.20\ \textrm {kg}\ \textrm {m}^{-3}$ for air and 
$\rho _{v} = 0.75\ \textrm {kg}\ \textrm {m}^{-3}$ for vapour. The water is modelled using the stiffened-gas equation of state
\begin{equation} p_{i} = \rho_{i}(\kappa_{i}-1)(e_{i}-e_{i,ref}) -\kappa_{i}{\rm \pi}_{\infty}, \end{equation}
where 
${\rm \pi} _{\infty }$ is yet another model parameter, and 
$\kappa _{w} = 2.35$, 
${\rm \pi} _{\infty } = 1$ GPa, 
$e_{w,ref} = 0\ \textrm {J} \textrm {kg}^{-1}$ and 
$\rho _{w} = 998\ \textrm {kg} \textrm {m}^{-3}$ (Le Métayer et al. Reference Le Métayer, Massoni and Saurel2004).
The interfacial pressure is defined as
\begin{equation} p_{I} =
\frac{\displaystyle \sum_{i=1}^{N} \left( p_{i}\sum_{j {\neq} i}^{N} \rho_{j}c_{j} \right)}{\displaystyle (N-1)\sum_{i =
1}^{N} \rho_{i}c_{i}},
\end{equation}
where 
$c_{i}$ is the speed of sound corresponding to the 
$i$th phase.
The pressure disequilibrium requires the internal-energy equation for each phase to be considered instead of the total energy equation of the mixture, while the conservation of the total energy of the mixture can still be written as
\begin{equation} \frac{\partial \rho E}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} [( \rho E + p ) \boldsymbol{u} + \boldsymbol{\varOmega} \boldsymbol{\cdot} \boldsymbol{u}] = 0 . \end{equation}Despite being redundant when computing the internal energy equation, the mixture total energy is also computed in practice to ensure numerical conservation of the total energy and consequently a correct treatment of shock waves (Saurel, Petitpas & Berry Reference Saurel, Petitpas and Berry2009).
2.2.2. Numerical procedure
The governing equations (2.1) are numerically solved using a splitting procedure, where the flow-associated terms on the left-hand side are dealt with separately from the terms associated with the relaxation procedure found on the right-hand side. The flow-associated terms are also split to first compute the hyperbolic terms using an explicit finite-volume Godunov scheme. The associated Riemann problem is computed using the Harten–Lax–van Leer–Contact (HLLC) approximate solver. Second, the surface tension terms are solved (Schmidmayer et al. Reference Schmidmayer, Petitpas, Daniel, Favrie and Gavrilyuk2017). Third, the right-hand side terms are solved using an infinite-relaxation procedure (Saurel et al. Reference Saurel, Petitpas and Berry2009) guaranteeing a unique pressure and providing a better estimate of the solution. Conservation of the total energy is ensured by correcting the non-conservative terms of the internal-energy equations through the total energy of the mixture (2.7).
Furthermore, a second-order-accurate MUSCL scheme is employed, where the time is integrated in two steps, the first acting as a predictor for the second (Schmidmayer et al. Reference Schmidmayer, Petitpas, Le Martelot and Daniel2020). The primitive variables are computed using a piecewise linear MUSCL reconstruction (Toro Reference Toro2013) with a monotonised central (Van Leer Reference Van Leer1977) slope limiter.
2.2.3. Computational problem, initial conditions and validation
In the experiment, sketched in figure 2(a), the laser-induced shock wave is initiated at a distance of approximately 5–7 mm from the stream of GEDs to enable sufficient shock wave energies to initiate the drop release. The laser-induced plasma also generates a vapour bubble that expands in the wake of the shock wave and can grow to a maximum radius of 2.5 mm, inducing a radial flow and pressure field that can influence the GED dynamics. The GEDs have a radius of approximately 1.2 mm and are separated from each other within the stream by a distance of approximately 4–6 mm, constrained by their generation method. The relatively short interbubble distances imply non-negligible effects on the individual GED dynamics and, to a weaker extent, shock reflections. At a later stage, the collapse of the vapour bubble induces a second shock wave, which further influences the drop mixing dynamics in the surrounding liquid. Considering the very short time scale with respect to the entire GED release dynamics, the laser-induced shock wave, the vapour bubble growth and collapse, and the second shock wave can technically be referred to as an effective single impulse.
Figure 2. (a) Schematic of the experimental conditions shortly after the shock wave inception. The region of interest (R.O.I.) indicates the region modelled numerically. An image sequence showing an example of shock–GED interaction is displayed with experimental radiographs in the upper halves and the corresponding ECOGEN simulation in the lower halves. The vapour and air bubbles are qualitatively displayed in blue and the shock waves in black. The indicated dimensionless time corresponds to that in panels (b–d). (b) Experimental, simulated and theoretical time evolution of the vapour bubble radius. The radius and time are respectively normalised to the theoretical maximum radius, 
$r_{vn,{max}}$, and the Rayleigh collapse time, 
$\tau _{c}$, given by the Keller–Miksis equation (2.8). Theoretical (c) pressure (2.9) and (d) velocity (2.10) fields induced by the vapour bubble motion at different distances. The pressure and velocity are normalised to the ambient pressure and characteristic speed of the flow, respectively. The thin lines delimit the range of normalised experimental radial distances from the vapour bubble centre, 
$r/r_{vn,{max}}$, and the bold lines show the mean.
The numerical framework is designed to simulate the experimental recordings of the dynamics of single GEDs interacting with laser-induced shock waves as faithfully as possible. Therefore, the air bubble, drop, laser-induced shock wave and the laser-induced vapour bubble are all modelled using three phases and their respective equations of state given in § 2.2.1. Water is applied to the entire domain and the drop, air is applied to the bubble and vapour to a small, highly pressurised nucleus, to represent the product of laser-induced plasma. The vapour nucleus induces the spherically propagating shock wave and grows similarly to a laser-induced vapour bubble. Only a single GED is simulated to allow for axisymmetry and mitigate the computational cost. Therefore, the air bubble, drop and vapour bubble centres are always located on the same axis to which a symmetrical boundary condition is applied, as depicted in figure 2(a). Non-reflecting boundary conditions are applied to the remaining sides of the numerical domain, which is made sufficiently large to prevent spurious waves induced at the non-reflecting boundaries from reaching the GED. Mesh stretching is applied outside of the region of interest (R.O.I.), where the air bubble, drop and vapour nucleus are located. A colour is assigned to the drop to enable tracking and thereby to reveal its release and dissemination into the surrounding liquid. The convergence of the bubble and drop dynamics is ensured by using a minimum of 104 grid points per drop diameter as shown in the grid dependence study of Appendix A. Solutions are saved at a frequency of 0.568 Mfps to match the experimental recording speed and to enable a direct comparison between the experiment and the numerical solution.
The shape and amplitude of the shock wave as well as the vapour bubble dynamics are solely driven by two initial numerical parameters, which are the initial radius, 
$r_{vn,0}$, and pressure, 
$p_{vn,0}$, of the vapour nucleus. An adequate combination of initial parameters 
$(r_{vn,0},p_{vn,0})$ is selected by relying on both the experimental hydrophone recordings of the shock wave pressure waveform and the high-speed recordings of the vapour bubble dynamics. The pair of initial parameters is estimated by iterating between the Keller–Miksis equation (Keller & Miksis Reference Keller and Miksis1980) and the power law describing the nonlinear pressure dissipation caused by the spherical spreading of the shock wave in the liquid (Vogel et al. Reference Vogel, Busch and Parlitz1996; Bokman et al. Reference Bokman, Biasiori-Poulanges, Meyer and Supponen2023b). The Keller–Miksis equation is employed to predict the combination of initial parameters 
$(r_{vn,0},p_{vn,0})$ that result in a good agreement between the simulated and experimental radius-time-curves of the vapour bubble. The equation reads
\begin{equation} \left( 1 - \frac{\dot{r}}{c_{0}}\right)r\ddot{r} + \frac{3}{2}\left( 1 - \frac{\dot{r}}{3c_{0}}\right)\dot{r}^{2} = \left( 1 + \frac{\dot{r}}{c_{0}}\right)\frac{ p_{b} - p_{0}}{\rho} + \frac{r}{\rho c_{0}}\frac{\mathrm{d}p_{b}}{\mathrm{d}t}. \end{equation}
Here, the bubble radius, interfacial velocity and interfacial acceleration are represented by 
$r$, 
$\dot {r}$ and 
$\ddot {r}$, respectively. The pressure on the liquid side of the bubble's interface is 
$p_{b} = p_{v} + p_{g,0}(r_{0}/r)^{3\kappa } - 2\gamma /r - 4\mu \dot {r}/r$, where the saturated vapour pressure is 
$p_{v}$ and an adiabatic formulation is used to model the non-condensible gases characterised by their partial pressure 
$p_{g,0}$, the bubble initial radius 
$r_{0}$ and the polytropic exponent of the vapour, 
$\kappa$. In (2.8), 
$(r_{0},p_{g,0})$ correspond to the driving pair parameter 
$(r_{vn,0},p_{vn,0})$. The surface tension, dynamic viscosity and speed of sound of the liquid are given by 
$\gamma$, 
$\mu$ and 
$c_{0}$ for water, respectively. However, the viscosity is set to zero in accordance with the numerical simulations that neglect it. For each combination of parameters which yield a good agreement between the Keller–Miksis equation and the experimental growth of the vapour bubble, the theoretical peak pressure at the location of the GED is verified by leveraging the power law describing the nonlinear pressure dissipation of the shock wave decaying as 
$p_{max}(r_{S})=p_{vn,0}( r_{vn,0}/r_{S})^{1.4}$ (Bokman et al. Reference Bokman, Biasiori-Poulanges, Meyer and Supponen2023b), where 
$r_{S}$ is the radius of the shock wave, and 
$r_{S}=0$ is the shock wave origin. The initial conditions agreeing well with the experimental shock wave and vapour bubble dynamics are 
$r_{vn,0}=25\ \mathrm {\mu }$m and 
$p_{vn,0} = 170$ GPa for the Keller–Miksis equation and 
$r_{vn,0}=25\ \mathrm {\mu }$m and 
$p_{vn,0} = 140$ GPa for the ECOGEN simulation, the results of which are discussed in Appendix B. Figure 2(b) displays a good agreement between the experimental recording, the Keller–Miksis equation and ECOGEN simulation up to 
$t/\tau _{c} = 0.45$, where 
$\tau _{c}=224\ \mathrm {\mu }$s is the Rayleigh collapse time corresponding to a maximum vapour bubble radius of 
$r_{vn,{max}}=2.44$ mm. The deviation beyond this instant is caused by the presence of a single GED (simulation) or multiple GEDs (experiment) acting as neighbouring free surfaces which shorten the oscillation time of the vapour bubble and causes its centre to move away from the free surfaces during the collapse phase, while producing a liquid jet in the same direction (Blake & Gibson Reference Blake and Gibson1981). The presence of multiple GEDs acting as free surfaces in the experiment has a greater impact on the vapour bubble dynamics as observed in the image sequence of figure 2(a) in comparison with the single GED present in the numerical simulation. The different vapour bubble dynamics in the experiment and simulation causes slight deviations in the GED dynamics. However, the numerical method is considered robust and a powerful tool capable of yielding additional insights on the behaviour of GEDs.
The pressure and velocity fields induced by the motion of the laser-induced vapour bubble affect the GED dynamics. For a spherical vapour bubble oscillation in a free field (Rayleigh Reference Rayleigh1917), the pressure in the liquid phase reads
\begin{equation} p_{vn}(r) = p_{0} + (p_{0}-p_{v})\left[ \frac{r_{vn}}{3r} \left( \frac{r_{vn, max}^{3}}{r_{vn}^{3}} -4\right) - \frac{r_{vn}^{4}}{3r^{4}} \left( \frac{r_{vn, max}^{3}}{r_{vn}^{3}} -1\right) \right], \end{equation}
at a distance, 
$r$, from the bubble centre. Note that at early times [
$t/\tau _{c} \rightarrow {O}(10^{-2})$], this simple model becomes significantly inaccurate. This is because (2.9) does not account for compressible effects related to the shock wave emission and neglects viscosity and the presence of non-condensible gases. In addition, the radial flow caused by the bubble motion in the surrounding liquid is
\begin{equation} u_{vn}(r) = \frac{r_{vn}^{2}\dot{r}_{vn}}{r^{2}}. \end{equation}
Here, 
$r_{vn}$ and 
$\dot {r}_{vn}$ are the vapour bubble radius and radial velocity as they evolve in time. These variables can be obtained from the Keller–Miksis solution displayed in figure 2(b) and used to give an estimate of the local pressure and radial flow velocity in the liquid surrounding the vapour bubble, plotted in figures 2(c) and 2(d), respectively, at the location, 
$r$, across the range covered by the experiments. The pressure is normalised to the ambient pressure and the velocity to the characteristic speed of the flow, 
$\sqrt {p_0/\rho }$. The pressure field surrounding the GEDs at their mean location from the vapour bubble centre, 
$r=5.6$ mm, is estimated to reach values 42 % lower than the ambient pressure once the vapour bubble reaches its maximum size resulting in an over-expansion of the air bubbles compared with solely shock-wave-driven dynamics. The corresponding radial flow velocity (maximum, 
$u_{vn}=1.03\ \textrm {m} \textrm {s}^{-1}$), on the other hand, may contribute to slightly displacing and deforming the GED. However, shock-wave-driven dynamics are expected to largely dominate the drop release process, their characteristic speeds being in the range 10–
$20~\textrm {m}~\textrm {s}^{-1}$.
3. Results and discussion
The observed dynamics resulting from the shock–GED interactions are generally similar to the numerical predictions of Biasiori-Poulanges et al. (Reference Biasiori-Poulanges, Bokman, Baumann and Supponen2022). Once the shock wave impacts the outer interface of the GED (at 
$t=0$), the latter effectively senses a strong pressure difference which initiates the bubble collapse. The time delay caused by the passage of the shock wave over the bubble span makes the collapse non-symmetrical, resulting in the proximal side to contract sooner and faster than the distal side. The bubble dynamics of the GED does not yield any noticeable differences to the dynamics of a simple shock–bubble interaction (Philipp et al. Reference Philipp, Delius, Scheffczyk, Vogel and Lauterborn1993; Ohl & Ikink Reference Ohl and Ikink2003; Johnsen & Colonius Reference Johnsen and Colonius2009; Biasiori-Poulanges et al. Reference Biasiori-Poulanges, Bokman, Baumann and Supponen2022) until the outer interface reaches the inner drop.
Depending on the relative location of the drop within the bubble, differences in the drop release processes are observed. They are herewith classified into three distinct regimes: drop-impact-, partial-deposition- and jet-impact-driven drop release.
3.1. The drop impact and partial deposition regimes
The drop impact regime yields weaker yet similar drop dynamics as those reported in the work of Biasiori-Poulanges et al. (Reference Biasiori-Poulanges, Bokman, Baumann and Supponen2022), and can be observed in the test case displayed in figure 3, which shows both numerical simulations and experimental radiographs. The initial bubble radius is 
$r_{b,0}=1.20$ mm and the drop-to-bubble size ratio is 
$\varepsilon =r_{d,0}/r_{b,0}=0.37$. The bubble centre is located 5.34 mm from the shock wave source and the drop rests near the proximal bubble side, 0.76 mm from the air shell's centre. The shock wave parameters upon contact with the GED are 
$p_{max} = 62.45$ MPa, 
$j = 2.79$ Pa s and 
$M=1.06$. Excellent agreement between the numerics and the experiment is observed if not for a slight over-expansion of the bubble at later stages in the numerical simulation.
Figure 3. Image sequence of an example of shock–GED interaction within the drop impact regime (
$We_{d} = 1940,\ Re_{d} = 11\ 000$). The upper halves display the numerically computed volume fractions of air and drop water in blue and black, respectively. The lower halves show corresponding experimental radiographs. The dimensionless time of the bubble, 
$t/(r_{b,0}\sqrt {\rho / p_{b,0}})$, starting from the shock contact with the GED, is indicated in the bottom-left corner. The inset figure shows the normalised bubble pressure map 
$p_{b}/p_{0}$.
The shock-induced collapse of the air bubble can be observed in frames 1–4 of figure 3 and is followed directly by a rebounding expansion phase in the subsequent frames. The drop's location on the bubble axis, very close to the proximal interface of the bubble, results in an instantaneous impact during the bubble contraction. A necking region then starts radially expanding from the contact point between the drop and the surrounding liquid, driving the latter into the air bubble in the form of a radially and axially propagating sheet jet, as can be observed on the second and third frames of figure 3. Such a sheet jet has been attributed to the strong water-hammer shock acting on the bubble cusp (Lesser & Field Reference Lesser and Field1983; Biasiori-Poulanges et al. Reference Biasiori-Poulanges, Bokman, Baumann and Supponen2022). However, the outer liquid impacts the drop displayed in figure 3 at an experimentally recorded impact speed of 
$u_{i,exp} = 12.62\ \textrm {m}\ \textrm {s}^{-1}$, close to the value of 
$u_{i,sim} = 10.98\ \textrm {m} \textrm {s}^{-1}$ retrieved from the numerical simulation along the axis shared between the bubble and drop. Here, the impact speed is defined as the relative velocity of the bubble interface, where the drop speed is typically negligible, being one to two orders of magnitude smaller. At such speeds, the drop is impacted at a Mach number in air of 
$M=0.04$, much lower than the limit of 
$M\geq 1$ at which compressibility effects should become relevant (Cheng, Sun & Gordillo Reference Cheng, Sun and Gordillo2022), resulting in a possible different driving mechanism for the sheet jet dynamics than the water hammer.
The splashing mechanism observed in figure 3 displays striking similarities with Worthington's first instantaneous photographs of a drop impacting a liquid surface (Worthington Reference Worthington1908) taken from A study of splashes in the early stages of necking, sheet jet inception and beyond. The Weber number, 
$We = 2\rho r_{d,0} u_{i}^{2} / \gamma$, describing the inertial to surface tension force balance, and the Reynolds number, 
$Re = 2\rho r_{d,0} u_{i} / \mu$, denoting the inertial to viscous force balance, are 
$We_{d} = 1940$ and 
$Re_{d} = 11\ 000$ for the drop in the experiment. The splash of a liquid drop on a flat liquid surface at such high Weber and Reynolds numbers adopts the shape of a crown as the sheet jet grows (Engel Reference Engel1966; Prosperetti & Oguz Reference Prosperetti and Oguz1993; Rein Reference Rein1996; Wang & Chen Reference Wang and Chen2000; Thoroddsen Reference Thoroddsen2002; Yarin Reference Yarin2006; Wang et al. Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023). This shape cannot be resolved herein; however, its consequences can be observed starting from the fourth frame of figure 3. As the thickened rim of the crown becomes unstable to the Rayleigh–Plateau instability, it forms a spray of drops which cannot be resolved individually in the radiographs, but can be observed upon clustering between the sheet jet and the bubble surface as shades of grey. In the simulation, the spray of small drops ejected from the crown tip can be seen stretching radially until it reaches the bubble interface, as highlighted on the fourth frame of figure 3. At the same moment, the sheet jet tip is redirected upstream after having impacted the bubble surface. However, due to the lower resolution of the sheet jet tip at the impact point, redirection of part of the jet downstream cannot be excluded. This has been observed in the splashing of drops on liquid surfaces Thoroddsen et al. (Reference Thoroddsen, Thoraval, Takehara and Etoh2011) and in the sheet jet cascade mechanism reported in violent shock-induced GED collapses (Biasiori-Poulanges et al. Reference Biasiori-Poulanges, Bokman, Baumann and Supponen2022). The sheet jet motion creates a distinctive region visible in the top corner of the fifth frame of figure 3. The inset displays the normalised pressure field, 
$p_{b}/p_{0}$, within the bubble, where a pressure drop of approximately 15 % is caused by fluid rotation. This rotating region likely contributes to the ejection of drops from the crown rim by vortex shedding (Jerome et al. Reference Jerome, Marty, Matas, Zaleski and Hoepffner2013). After ejection, the impact of these drops against the bubble wall can be observed in frames 5 and 6 where they disrupt the initially smooth air–water interface. This is especially visible in the highlighted region of the experimental part of the sixth frame of figure 3. Note that slight surface instabilities are also visible in all numerical simulations. However they are found to be grid dependent and are considered a numerical artifact with no incidence on the present work. Note that the bubble cusp located at the base of the sheet jet is visible in the radiographs as a vertical line starting from the fifth frame of figure 3, where the x-rays integrate all phase discontinuities along the beam path. It is, however, not visible in the two-dimensional representation of the axisymmetric numerical simulation.
In frames 7–9, the rim of the crown proceeds towards the central axis of the GED, and ends up separating the bubble into two large volumes of air. A central ‘twin’ jet forms, travelling both upstream and downstream and eventually piercing both sides of the bubble as revealed by the simulations. Although part of the drop is entrained by the sheet jet, its main core barely moves as it gets released locally during the overall process of bubble compression. The air is simply transferred to the distal side of the drop as it coalesces with the surrounding liquid. The central jet travelling upstream impacts the drop once more, increasing its deformation and dissemination. Note that a small fraction of drop contents is entrained by the sheet jet, which suggests that the splashing process contributes to the dissemination of the drop.
The partial deposition regime is shown in figure 4, where the numerical simulation, indicating the volume fraction of air in grey and drop in black, is displayed in the upper half of the image and the corresponding experimental radiograph is displayed in the lower half. Here, the highly deformed initial shape of the bubble in the experiment makes a direct comparison of the entire gas-encapsulated drop irrelevant. The liquid jet (not shown here) forms at the point of highest interfacial curvature in the experiment, which differs from the axisymmetric simulation results, where the bubble is a sphere. The full experimental and simulated GED dynamics are available in the supplementary movies are available at https://doi.org/10.1017/jfm.2024.1124. Since the drop dynamics take place at the distal side of the bubble, figure 4 solely focuses on that part.
Figure 4. Image sequence of an example of shock–GED interaction within the partial deposition regime (
$We_{d} = 16,\ Re_{d} = 777$). The upper halves display the numerically computed volume fractions of air and drop water in grey and black, respectively. The lower halves show corresponding experimental radiographs. The dimensionless time, 
$t/(r_{b,0}\sqrt {\rho / p_{b,0}})$, starting from the shock contact with the GED, is indicated in the bottom-left corner.
The initial size of the bubble is 
$r_{b,0}=1.14$ mm and the drop-to-bubble size ratio is 
$\varepsilon =0.24$. The bubble centre is located 5.67 mm from the shock wave origin and the drop centre 0.82 mm downstream from the bubble centre. The shock wave parameters are 
$p_{max}=55.93$ MPa, 
$j=2.50$ Pa s and 
$M=1.06$. The drop contacts the outer liquid at a relatively low experimentally recorded speed of 
$u_{i, exp}=1.45\ \textrm {m}\ \textrm {s}^{-1}$, close to that extracted from the numerical simulation of 
$u_{i, sim}=1.65\ \textrm {m}\ \textrm {s}^{-1}$. Since the drop is located near the distal side of the bubble, the closest interface is the slowest part of the bubble during its asymmetrical collapse. The Weber and Reynolds numbers associated with the drop upon impact in the experiment are 
$We_{d}=16$ and 
$Re_{d}=777$.
Here, the physical problem is otherwise similar to that presented in the drop impact regime, except for the weaker inertia as suggested by the dimensionless numbers which are two orders of magnitude lower. In particular, the associated Weber number suggests capillary effects, while viscosity is still expected to play a minor role due to the large Reynolds number. The drop dynamics observed in figure 4 display strong similarities with the deposition and partial coalescence of drops, reported in an abundant number of works (Charles & Mason Reference Charles and Mason1960; Blanchette & Bigioni Reference Blanchette and Bigioni2009; Kavehpour Reference Kavehpour2015; Kulkarni et al. Reference Kulkarni, Lolla, Tamvada, Shirdade and Anand2021). However, although the drop pinch-off process is governed by surface tension, the underlying mechanism here appears fundamentally different. Indeed, the partial coalescence of drops impacting flat liquid surfaces at very low speeds and the resulting ejection of daughter droplets have been attributed to the motion of capillary waves on the drop surface, generated by the formation of a neck at the drop–liquid interface in the early stages of coalescence (Blanchette & Bigioni Reference Blanchette and Bigioni2006). The momentum of sufficiently strong waves has been reported to distort the drop significantly when converging towards its summit, stretching the drop away from the liquid surface and inducing a collapse in the direction perpendicular to the stretching, which results in the production of a daughter droplet. Here, the deformation of the drop does not appear to be driven by capillary waves, but rather by the motion of the bubble interface following its interaction with the shock wave.
Once the bubble collapse initiates, the air film between the drop and the outer liquid becomes thinner. As soon as both liquids contact, the deformation of the drop and the coalescence begin. Similarly to the drop impact regime, a radially expanding necking region can be observed at the liquid–liquid interface, as highlighted in the second frame of figure 4. The early stage of the drop dynamics follows the ones reported during partial coalescence processes up to the end of the collapse of the air bubble surrounding the drop. Once the bubble rebounds, the reversed motion of its interface coupled with the capillary force linking the drop to the outer liquid deforms the drop into a pear-like shape, as observed in the third frame. Unlike in previously reported partial coalescence processes, capillary waves are absent and the drop deformation solely takes place at the proximal side of the drop, while its distal side remains at its initial location. The oscillatory motion of the bubble interface, as it pulls on the drop during the bubble expansion, drives the collapse of the newly formed neck by stretching the drop, as observed in the two last frames of figure 4. In the current work, however, daughter droplets rarely fully detach due to the second shock wave caused by the collapse of the laser-induced vapour bubble driving a second collapse of the air bubble. This often happens before a daughter droplet has formed completely making the drop dynamics difficult to follow after that point.
The early stage of the drop impact and partial deposition regimes are, in fact, different inertial regimes of drops being impacted by a liquid surface, which is fairly well understood and can be described through scaling laws. The spreading radius of the necking region between the drop and the outer liquid has previously been reported to scale with time as 
$r_{n}/r_{d,0} \propto \sqrt {2u_{i}t/r_{d,0}}$ (Josserand & Zaleski Reference Josserand and Zaleski2003), where 
$r_{d,0}/u_{i}$ is the characteristic impact time of the drop, based on the truncated sphere approximation and solely on geometric considerations. However, such approximations do not respect the continuity equation. Riboux & Gordillo (Reference Riboux and Gordillo2014) instead found using Wagner's theory (Wagner Reference Wagner1932),
\begin{equation} \frac{r_{n}}{r_{d,0}}=\sqrt{\frac{3u_{i}t}{r_{d,0}}}, \end{equation}
to be a more accurate formulation, while still neglecting viscosity or air entrapment. The model is valid until the necking radius reaches the drop's initial size, 
$r_{n}/r_{d,0} \leq 1$, although good agreement has also been reported at later times (Josserand & Zaleski Reference Josserand and Zaleski2003; Howland et al. Reference Howland, Antkowiak, Castrejón-Pita, Howison, Oliver, Style and Castrejón-Pita2016). Figure 5 shows the temporal evolution of the spreading radius of the necking region recorded in the experiment and simulation of the drop impact and partial deposition regimes against the model given by (3.1). Excellent agreement is found between the numerical simulations and the model for the entire range of validity of the drop impact regime. In the partial deposition regime, however, the neck radius diverges from the model at 
$u_{i}t/r_{d,0} = 0.18$, decreasing as the air bubble expands and stretches the drop, as observed in the three last frames of figure 4. The experiment follows the same trends as the numerical simulations, consistent with the excellent agreement reported in figures 3 and 4. The slightly higher values are likely caused by the drop being deformed and already pressed against the bubble interface at the moment of the impact in the drop impact regime, while the simulation is initiated with both the drop and bubble being modelled as spheres. They are, however, well within the experimental uncertainty. A fitted version of (3.1), 
$r_{n}/r_{d,0}=C\sqrt {3u_{i}(t-t_{\mu })/r_{d,0}}$, has been suggested (Jian et al. Reference Jian, Josserand, Popinet, Ray and Zaleski2018; Wang et al. Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023), where 
$C$ and 
$t_{\mu }$ are fitting parameters representing a scaling constant and the time delay correcting for viscous and cushioning effects, respectively. The fitted parameters found for the experiment and simulation are 
$C_{exp}=1.03$, 
$C_{sim}=1.02$, 
$t_{\mu,exp}=-1.78\ \mathrm {\mu }$s and 
$t_{\mu,sim}=0\ \mathrm {\mu }$s for the drop impact regime and 
$C_{exp}=1.11$, 
$C_{sim}=1.06$ and 
$t_{\mu,exp}=-1.83\ \mathrm {\mu }$s, 
$t_{\mu,sim}=0.53\ \mathrm {\mu }$s for the partial deposition regime. The scaling constants are for all cases close to one, which supports the modelling of the necking region. The very short time delays suggest that air cushioning and viscous effects are negligible. The negative values are likely caused by the necking region having begun forming at an earlier time in the experiment, but the exact timing of impact cannot be resolved at the current spatiotemporal resolution. The fitted curves are not displayed in figure 5 to maintain visual clarity since they are so close to the uncorrected model.
Figure 5. Temporal evolution of the necking radius for the drop impact (
$We_{d} = 1940$) and partial deposition (
$We_{d} = 16$) regimes. All quantities are normalised to the drop's initial radius, 
$r_{d,0}$, the drop impact speed, 
$u_{i}$ and the drop characteristic time, 
$r_{d,0}/u_{i}$.
In the drop impact regime, the necking region eventually evolves into a sheet jet starting from the second frame of figure 3, which corresponds to 
$u_{i}t/r_{d,0}=0.51$. Since the sheet jet entrains a small fraction of drop contents with it, understanding and predicting the way it spreads within the air bubble is of interest. The model given by (3.1) agrees well with the temporal evolution of the spreading radius of the sheet jet beyond the limit of 
$r_{n}/r_{d,0} \leq 1$, up to 
$u_{i}t/r_{d,0}=0.69$ as displayed in figure 6(a), where 
$r_{s}=r_{n}$ before the apparition of the sheet jet. Past that time, the spreading radius of the sheet jet converges to a more or less constant value of 
$r_{s}/r_{d,0} = 1.5$ up to 
$u_{i}t/r_{d,0}=2.75$, corresponding to the second row of images in figure 3. This is likely caused by the sheet jet spraying daughter droplets from its tip as it becomes unstable starting from 
$u_{i}t/r_{d,0} \approx 1$. Note that a whole-process theory for drop impact, such as that proposed by Bisighini et al. (Reference Bisighini, Cossali, Tropea and Roisman2010) for the spreading of the drop crater, does not capture the spreading of the sheet jet because of the shock-induced bubble oscillation. Indeed, during the bubble collapse and expansion, Bisighini et al.'s (Reference Bisighini, Cossali, Tropea and Roisman2010) theory over- and under-estimates the sheet jet spreading, respectively. The radial speed of the neck and early sheet jet can be estimated by differentiating (3.1) with respect to time, yielding
\begin{equation} \frac{\dot{r}_{n}}{u_{i}}=\frac{3}{2\sqrt{3u_{i}t/r_{d,0}}}. \end{equation}
This provides a good approximation of the spreading speed of the neck at an early stage as displayed in figure 6(b). However, the model soon overestimates the recorded radial speed as the radial evolution of the sheet jet halts. The radial deceleration occurs prior to the onset of instabilities at the sheet jet tip which eject daughter droplets and is likely the main cause of them. The maximum speed is observed during the expansion of the necking region, with 
$\dot {r}_{s} = 28.02\ \textrm {m}\ \textrm {s}^{-1}$. Past 
$u_{i}t/r_{d,0} = 2.75$, the sheet jet radius starts decreasing and the radial speed becomes negative as the sheet jet tip converges back towards the central axis of the GED, as observed in the last row of frames of figure 3.
Figure 6. Temporal evolution of the crown's (a) radius, (b) radial speed, (c) height and (d) speed along the GED axis (
$We_{d} = 1940, Re_{d} = 11\ 000$). All quantities are normalised to the drop's initial radius, 
$r_{d,0}$, the drop impact speed, 
$u_{i}$ and the drop characteristic time, 
$r_{d,0}/u_{i}$.
The height of the splash, 
$h_{s}$, measured from the drop proximal side, and speed at which it travels downstream are measured and displayed in figures 6(c) and 6(d), respectively. Prior to trackable sheet jet formation, 
$h_{s}$ is measured as the distance of the cusp of the bubble to the initial position of the drop's proximal side. One notices a slight acceleration phase during the necking stage up to sheet jet formation, where the sheet jet reaches a maximum speed of 
$\dot {h}_{s}=33.44\ \textrm {m}\ \textrm {s}^{-1}$ at 
$u_{i}t/r_{d,0} = 0.36$. The splash then grows at a more or less constant speed of 
$\dot {h}_{s,exp}\approx 15.72\ \textrm {m}\ \textrm {s}^{-1}$ and 
$\dot {h}_{s,sim}\approx 15.31\ \textrm {m}\ \textrm {s}^{-1}$ until it reaches the distal side of the bubble at 
$u_{i}t/r_{d,0} \approx 2.8$. Quantitative analysis beyond that point becomes arduous in the experiment due to the strong deformations of the bubble interface.
3.2. The jet impact regime
The third regime of shock-induced GED release, named the jet impact regime and reported for the first time herein, is shown in figure 7. Despite the considerable initial deformation of the bubble at the time of shock wave impact in the experiment, the numerical simulation agrees reasonably well with the experiment.
Figure 7. Image sequence of an example of shock–GED interaction within the jet impact regime (
$We_{d} = 2810$, 
$We_{j} = 1580$, 
$Re_{d} = 9730$, 
$Re_{j} = 5450$). The upper halves display the numerically computed volume fractions of air and drop water in red and black, respectively. The lower halves show corresponding experimental radiographs. The dimensionless time, 
$t/(r_{b,0}\sqrt {\rho / p_{b,0}})$, starting from the shock contact with the GED, is indicated in the bottom-left corner.
In this regime, the properties of the shock wave and the GED are very similar to those presented in the previous regimes. The initial bubble radius is 
$r_{b,0}=1.20$ mm and the drop-to-bubble size ratio is 
$\varepsilon =0.19$. The bubble centre is located 5.83 mm from the shock wave source. The shock wave parameters are 
$p_{max}=54.45$ MPa, 
$j=2.44$ Pa s and 
$M = 1.06$. The shock-induced bubble dynamics in frames 1–3 are unaffected by the small drop located 0.26 mm downstream from the centre of the bubble, and thus avoiding any early contact with the bubble wall. As a result, the proximal side of the bubble has sufficient space to evolve into a liquid jet which can be observed in the following frames. Such high-speed jets are expected after the interaction of bubbles with shock waves of substantial pressure impulse (Philipp et al. Reference Philipp, Delius, Scheffczyk, Vogel and Lauterborn1993; Ohl & Ikink Reference Ohl and Ikink2003; Bokman et al. Reference Bokman, Biasiori-Poulanges, Meyer and Supponen2023b). Here, jetting is likely favoured by weak surface tension effects associated with the large size of the bubbles, as well as the additional contribution of the flow induced by the growth of the laser-induced vapour bubble. This might explain the discrepancies between the observed and computed speeds of the jet since the vapour bubble dynamics is slightly different in the experiment and simulation as displayed in figure 2(b). In the fourth radiograph, the jet impacts the liquid drop at a speed of 
$u_{i, exp}=21.02\ \textrm {m} \textrm {s}^{-1}$, which is faster than in the numerical simulation, 
$u_{i, sim}=16.96\ \textrm {m} \textrm {s}^{-1}$. To the best of the authors’ knowledge, only a few studies have investigated the impact of liquid jets with drops (Quetzeri-Santiago & Rivas Reference Quetzeri-Santiago and Rivas2023; Mulbah et al. Reference Mulbah, Kang, Song and Yin2024). Previous investigations on binary drops (Qian & Law Reference Qian and Law1997; Pan & Suga Reference Pan and Suga2005; Pan, Chou & Tseng Reference Pan, Chou and Tseng2009; Roisman et al. Reference Roisman, Planchette, Lorenceau and Brenn2012; Cimpeanu & Moore Reference Cimpeanu and Moore2018), however, show strong similarities with the dynamical process observed in the jet impact regime. One can simplify the liquid jet impact with a drop as an unequal-sized drop collision (Ganti, Khare & Bravo Reference Ganti, Khare and Bravo2020), where an almost-cylinder-shaped drop impacts a larger spherical liquid drop. The present analogy is restricted to initial jet radii smaller or equal to the radii of the drops they are impacting, 
$r_{j,0} \leq r_{d,0}$. The smaller size of the jet with respect to the drop induces drop deformations upon impact in the shape of a thin sheet jet opening to the opposite direction of the travelling jet, i.e. upstream. Here, the experimental Weber and Reynolds number computed for the drop and jet are 
$We_{d} = 2810$, 
$We_{j} = 1580$ and 
$Re_{d} = 9730$, 
$Re_{j} = 5450$, which are in the same order of magnitude as the inertia-driven impact in the drop impact regime. At the early stage, a necking region forms, starting from the fourth frame in figure 7, which evolves into a sheet jet that gives a bag-like shape to the drop. As the bag spreads within the air bubble, it adopts an inflection point near its rim in the numerical simulation, with the sheet jet tip pointing perpendicularly to the jet axis, while in the experiment, the rim of the bag points upstream, parallel to the jet axis. This is particularly visible in the sixth frame of figure 7. The higher impact velocity in the experiment certainly contributes to this discrepancy, however other factors such as the assumption of two-dimensional axisymmetry in the simulation might also play a role. Due to the high momentum of the jet compared with the small liquid mass of the drop, the jet continues to travel downstream, displacing and deforming the drop located on its path. At first, the drop mainly deforms while maintaining its distal side at the same location, with most of its mass being pushed to the sides of the jet. This can be observed as the spreading radius of the sheet jet expands and the drop thickness along the jet axis reduces, as seen in frames 4–6 of figure 7. Once the drop has adopted a fine bag-like structure with a large surface area, its distal side starts moving with the jet towards the back of the bubble. The drop moves faster in the experiment than in the numerical simulation, which can easily be explained by the greater experimental impact speed. On the two last frames of figure 7, the bag has become so thin that it is arduous to track both in the experiment and numerical simulation. Finally, the drop impacts the distal bubble surface and mixes with the surrounding liquid, as seen in the last frame of figure 7, provoking a slight increase in pressure at the impact point. The drop rim also experiences instabilities turning it into a finger-like structure from which some liquid is ejected in the form of a spray of daughter droplets. This is very similar to the crown instability reported in the drop impact regime.
The spreading of the necking region and early formation of the sheet jet is quantified in figure 8(a). Here quantities are non-dimensionalised to the initial radius of the jet tip at the moment of impact, 
$r_{j,0}$, instead of the initial drop radius. The scaling of quantities by the jet's initial radius, where 
$r_{j,0}= 0.13$ mm, is herein used because of the smaller jet size with respect to the drop, in which it becomes the impactor and the drop becomes the impacted. This new scaling is further justified by the simulated dynamics of drops of different sizes (
$\varepsilon =0.19$ and 
$\varepsilon =0.36$) resting at the same location in a 1.20-mm bubble approximately collapsing on a single curve. The initial jet radius and impact speed for the larger drop are 
$r_{j,0}= 0.15$ mm and 
$u_{i}= 14.03\ \textrm {m}\ \textrm {s}^{-1}$, respectively. Unsurprisingly, the spreading of the necking region displays once more a square root scaling to the dimensionless time. However, a better agreement through
\begin{equation} \frac{r_{n}}{r_{j,0}}=\sqrt{\frac{2u_{i}t}{r_{j,0}}}, \end{equation}
compared with the corrected factor used in (3.1), is shown in figure 8(a). This is likely caused by the cylindrical shape of the impacting jet and, in fact, a similar scaling has previously been reported for the splashing of a jet on a thin liquid film (Bokman et al. Reference Bokman, Biasiori-Poulanges, Lukić, Bourquard, Meyer, Rack and Supponen2023a). A best fit of 
$r_{n}/r_{j,0}=C\sqrt {2u_{i}(t-t_{\mu })/r_{j,0}}$, where 
$C$ and 
$t_{\mu }$ are fitting parameters, to the experimental case and both numerical simulations yield 
$C_{exp}=0.92$, 
$C_{sim(\varepsilon =0.19)}=0.97$, 
$C_{sim(\varepsilon =0.36)}=1.00$ and 
$t_{\mu }=0\ \mathrm {\mu }$s for all. The scaling constants are close to one for both simulations, and larger for the experiments, which is expected considering the large uncertainty at such small scales. The absence of a time delay suggests once more the absence of air cushioning and viscous effects.
Figure 8. Temporal evolution of (a) the necking radius, the drop's (b) opening radius, (c) radial speed, (d) location and (e) speed along the GED axis from the first instant the jet contacts the drop. All quantities are normalised to the jet's initial radius, 
$r_{j,0}$, the drop impact speed, 
$u_{i}$ and the jet and drop characteristic times, 
$r_{j,0}/u_{i}$ and 
$r_{d,0}/u_{i}$. The drop-to-bubble size ratio, 
$\varepsilon$, is indicated for a 1.20-mm bubble.
In the context of payload delivery, the deformation and displacement of the drop are relevant. The spreading radius of the drop rim, 
$r_{d}$, is substantially smaller in the experiment than in the simulations, as shown in figure 8(b), and is expected from the differences in shape observed in figure 7. A square root scaling with the dimensionless time is obvious and, interestingly, the experimental spreading radius agrees with the model given by (3.3), also shown in figure 8(b). The speed at which the drop rim opens, 
$\dot {r}_{d}$, can be estimated from the spreading of the neck, assuming 
$r_{d}=r_{n}$. Differentiating (3.3) with respect to time yields
\begin{equation} \frac{\dot{r}_{n}}{u_{i}}=\frac{1}{\sqrt{4u_{i}t/r_{j,0}}}. \end{equation}
The model agrees reasonably well with the experimentally observed radial speeds and slightly underestimates the numerical results, as observed in figure 8(c). At later stages, the agreement between the experimental speed and numerical speed improves for 
$\varepsilon =0.19$, where 
$\dot {r}_{d, exp} = 4.21\ \textrm {m}\ \textrm {s}^{-1}$ and 
$\dot {r}_{d, sim(\varepsilon =0.19)} = 4.24\ \textrm {m}\ \textrm {s}^{-1}$ at 
$u_{i}t/r_{j,0}\approx 19$. The radial speed for the larger drop at a similar dimensionless time is 
$\dot {r}_{d, sim(\varepsilon =0.36)} = 1.19\ \textrm {m}\ \textrm {s}^{-1}$.
The displacement of the drop along the jet axis is another parameter of interest in the context of payload delivery. Before being set in motion, the distal side of the liquid drop stays immobile for 
$u_{i}t/r_{d,0} \lesssim 2$, where 
$r_{d,0}/u_{i}$ is the characteristic time of the drop and 
$u_{i}t/r_{d,0} = 0$ is the time at which the jet contacts the drop. Here, Buckingham–Pi theory is employed to find an equation to describe the displacement of the distal point on the drop surface, 
$x_{d}$, once it is set in motion by the jet. The theorem yields three different groups based on the main quantities of interest, mainly 
${\rm \pi} _{1} = x_{d}/r_{d,0}$, 
${\rm \pi} _{2} = r_{j,0}/r_{d,0}$ and 
${\rm \pi} _{3} = u_{i}t/r_{d,0}$, which can be expressed as 
${\rm \pi} _{1} = f({\rm \pi} _{2},{\rm \pi} _{3})$ in the form of
\begin{equation} \frac{x_{d}}{r_{d,0}} =
D \left(\frac{r_{j,0}}{r_{d,0}}\right)^{\alpha_{1}}
\left[\frac{u_{i}(t-t_{x_{d}})}{r_{d,0}}
\right]^{\alpha_{2}}, \end{equation}
where 
$D$, 
$\alpha _{1}$ and 
$\alpha _{2}$ are scaling constants to be fitted. A best fit from the first instant the distal point is set in motion, 
$t_{x_{d}}$, corresponding to 
$t_{x_{d}, exp}=33.44\ \mathrm {\mu }$s in the experiment and 
$t_{x_{d}, sim(\varepsilon =0.19)}=36.96\ \mathrm {\mu }$s and 
$t_{x_{d}, sim(\varepsilon =0.36)}=73.92\ \mathrm {\mu }$s in the numerical simulations, until 
$u_{i}t/r_{d,0}=10$ yields 
$D_{exp} = 0.13$, 
$D_{sim(\varepsilon =0.19)} = 0.14$ and 
$D_{sim(\varepsilon =0.36)} = 0.12$ and 
$\alpha _{2} = 1.5$ for all data sets. This yield the following scaling law
\begin{equation} \frac{x_{d}}{r_{j,0}} = D \left[\frac{u_{i}(t-t_{x_{d}})}{r_{d,0}} \right]^{3/2}, \end{equation}which adequately predicts the displacement of the distal side of the drop as shown in figure 8(d). Based on (3.6), an estimation of the speed at which the drop comes into contact with the bubble distal side can be drawn by simple time differentiation:
\begin{equation} \frac{\dot{x}_{d}}{u_{i}} = D \left(\frac{r_{j,0}}{r_{d,0}}\right) \left[\frac{u_{i}(t-t_{x_{d}})}{r_{d,0}} \right]^{1/2}. \end{equation}
Figure 8(e) shows the time evolution of the speed of the drop's distal side. A relatively good agreement with the model in (3.7) using the fitted parameters 
$A$ is obtained for each case. The terminal speed at which the drop travels right before impacting the distal bubble interface is 
$\dot {x}_{d, exp} = 7.36\ \textrm {m} \textrm {s}^{-1}$ in the experiment, 
$\dot {x}_{d, sim(\varepsilon =0.19)} = 5.94\ \textrm {m}\ \textrm {s}^{-1}$ in the corresponding simulation and 
$\dot {x}_{d, sim(\varepsilon =0.36)} = 2.76\ \textrm {m}\ \textrm {s}^{-1}$ for the larger drop. Despite a non-negligible displacement of the part of the drop directly located on the jet's path towards the back of the bubble, the entire drop's centre of mass, computed from the numerical simulations, does not evolve significantly as observed in figure 8(d). This is mainly caused by the lateral displacement of the drop rim perpendicular to the jet counterbalancing the drop displacement towards the distal bubble side.
Note that some regimes can also be combined. For example, a drop located on the GED axis, close to the distal side might first coalesce and then be impacted by a high-speed liquid jet. This leads to a more complex type of drop release which is outside of the scope of the current work.
3.3. Off-centre collision
Only axisymmetric cases have thus far been investigated, as they allow for a direct comparison with the numerical simulations, which are kept axisymmetric to minimise computational cost. However, in practice, the initial location of a drop within the bubble is random at the moment the shock wave impacts the GED, and the interfacial collisions are often asymmetric and off-centre. In the drop impact and jet impact regimes, the dynamics are affected by off-centre collisions, while in the partial deposition regime their effect is minor. The drop's eccentricity with respect to the impact is quantified by the impact parameter, 
$\psi _{d} = y_{d,0}/r_{d,0}$, which is defined as the distance of the drop's centre to the central axis of the bubble normalised to the initial radius of the drop.
In the drop impact regime, off-centre collisions induce asymmetrical sheet jets, as observed in the example of figure 9(a). The impact parameter here is 
$\psi _{d} = 0.42$ and the initial drop-to-bubble size ratio is 
$\varepsilon =0.40$. The shock wave inducing the GED dynamics has a peak pressure, pressure impulse and Mach number of 
$p_{max}=53.38$ MPa, 
$j=2.39$ Pa s and 
$M=1.06$, respectively. The asymmetric sheet jet dynamics shows strong similarities with the single-sided splashing reported in oblique drop impacts on liquid surfaces (Gielen et al. Reference Gielen, Sleutel, Benschop, Riepen, Voronina, Visser, Lohse, Snoeijer, Versluis and Gelderblom2017). In the present study, the bubble collapse is asymmetrical, with the proximal side collapsing first and faster than its distal side. The part of the drop located on the bubble axis, where the bubble collapse speed is maximal, is impacted by the bubble wall at a higher velocity than the lower part of the drop, which results in a more violent splashing on the upper side of the drop. Despite the asymmetric sheet jets, in essence, the splashing dynamics stay fundamentally the same as observed in the axisymmetric drop impact regime.
Figure 9. Radiographs showing both the initial location of a drop in a bubble and the off-centre (a) drop impact (
$We_{d} = 768$, 
$Re_{d} = 7180$) and (b) jet impact (
$We_{d} = 631$, 
$We_{j} = 253$, 
$Re_{d} = 14\ 600$, 
$Re_{j} = 5860$) regimes. The dimensionless time, 
$t/(r_{b,0}\sqrt {\rho / p_{b,0}})$, starting from the first contact with the surrounding liquid, drop and jet impact parameter, 
$\psi _{d}$, and, 
$\psi _{j}$, and jet deflection angle, 
$\theta _{j}$, are indicated. The central axis of the bubble and drop centre are shown as a dotted line and dark circle, respectively. (c) Jet deflection angle evolution as a function of the impact parameter.
The off-centre impact dynamics of a drop located in the jet impact regime is shown in figure 9(b). The impact parameter here is slightly modified to account for the thickness of the impacting jet, and is defined as 
$\psi _{j} = y_{d,0}/(r_{j,0}+r_{d,0}) = 0.58$. The initial drop-to-bubble size ratio is 
$\varepsilon =0.41$ and the shock wave has a peak pressure, pressure impulse and Mach number of 
$p_{max}=57.36$ MPa, 
$j=2.57$ Pa s and 
$M=1.06$, respectively. The off-centre jet impact on the drop results in the jet being deflected from its original course. Here, the deflection angle of the jet, 
$\theta _{j}$, can be estimated solely based on the impact parameter, 
$\psi _{j}$, as displayed in figure 9(c). Results indicate that upon impact, the jet deflection angle is dictated by the local inclination of the drop surface at the contact point with respect to the central drop axis, and can be expressed as
\begin{equation} \theta_{j} = \arccos{ ( \psi_{j} )} = \arccos{ \left( \frac{y_{d,0}}{r_{d,0}+r_{j,0}} \right)}. \end{equation}
A good agreement between the model and the experimentally observed deflection angles is found. This suggests a certain degree of rigidity of the liquid drop upon contact. In fact, the drop acts as an unmovable rigid surface which reorients the jet completely in the first instants after the impact. At a later stage, drop deformation comes into play, which is why only the tip of the jet is deflected by 
$\theta _{j}$ and the rest of the drop adopts a more complex shape. The tip of the jet, which is likely mixed with drop contents, finally impacts the bubble interface and mixes with the surrounding liquid.
3.4. Drop dynamics regime map
The two most important parameters dictating the drop release process in GEDs are the initial location of the drop within the bubble and the shock wave energy that results in different collapse speeds. The different regimes are mapped in figure 10(a) based on the initial location of the drop within the bubble at the moment the shock wave arrives. The shock parameters are kept approximately constant with a mean shock wave peak pressure and pressure impulse of 
$p_{max} = 59.11 \pm 13.22$ MPa and 
$j = 2.65 \pm 0.6$ Pa s, respectively. The 
$x$–
$y$-coordinates of the initial location of the drop are normalised to the bubble's initial radius, 
$r_{b,0}$. The drop-to-bubble size ratio, 
$\varepsilon$, is also provided, showing that in most cases, the drop size is 40 % of the bubble size and is only smaller in some isolated cases. The consistent drop-to-bubble size ratio is a result of the high repeatability of the GED production.
Figure 10. (a) Mapping of the initial position of the drop within the bubble, normalised to the bubble's initial radius, 
$r_{b,0}$, as it interacts with a shock wave of pressure impulse 
$j = 2.65 \pm 0.6$ Pa s. The drop release regimes, namely the drop impact regime (blue), partial deposition regime (grey) and jet impact regime (red), are indicated for a range of drop-to-bubble size ratios, 
$\varepsilon = r_{d,0}/r_{b,0}$, which is qualitatively emphasised by the marker size. The dashed circle shows the minimum equivalent size reached by the bubble. (b) The local curvature of the bubble at the impact location, 
$\mathcal {K}_{b}$, normalised to the drop's curvature, 
$\mathcal {K}_{d}$, and the drop's Weber and Reynolds numbers, delimit a regime map for all three drop release regimes. Quantities are collected at the moment of impact. The vertical and horizontal limits are 
$We^{1/2}Re^{1/4} = 76$ and 
$\mathcal {K}_{b}/\mathcal {K}_{d} = 1$, respectively.
The release of a drop located on the proximal side of the bubble takes place within the drop impact regime. That is because the proximal side of the bubble undergoes higher accelerations and speeds compared with the rest of the bubble as it collapses asymmetrically. This region appears to span along the entire 
$y/r_{b,0}$ domain as long as 
$x/r_{b,0} \leq -0.1$, where 
$(x/r_{b,0},y/r_{b,0}) = (0,0)$ is the bubble centre. In this region, a liquid jet does not have time to form and the drops are impacted by a relatively flat or even concave surface. The partial deposition regime instead occupies the annular region in the range 0.6–1.0 of the bubble's initial radius for 
$x/r_{b,0} > -0.1$. In this region close to the bubble wall, the interface contracts less violently than the proximal bubble side, where the jet forms. The mean of all the minimum equivalent bubble radii, 
$r_{min}=0.79r_{b,0}$, recorded from the bubble's observable area in the experiment after their contraction, is shown as a dashed line in figure 10(a). The drops located in the region between this circle and the one denoting the bubble's initial radius inevitably contact the outer liquid, initiating the partial coalescence regime. Larger drops may be located closer to the bubble centre and still coalesce. This is highlighted by the higher number of large drops near the bubble centre in comparison with the smaller drops located closer to the bubble wall. The jet impact regime only occurs in the central region of the GED where the jet has space to form and hits the drop before coalescence with other parts of the bubble wall can take place. The drops in this regime are located relatively close to the central bubble axis. This region is within the range 
$x/r_{b,0} = -0.1$ to 
$x/r_{b,0} = 0.5$ and 
$y/r_{b,0} = -0.5$ to 
$y/r_{b,0} = 0.5$.
The identified regime map in figure 10(a) is, however, rather specific to the shock type and bubble sizes used here. Below a certain shock pressure threshold, bubbles only respond through radial oscillations without producing any jets (Ohl & Ikink Reference Ohl and Ikink2003; Bokman et al. Reference Bokman, Biasiori-Poulanges, Meyer and Supponen2023b). Under such conditions, the regime map would look different, as the jet impact regime would be absent from the map. Similarly, smaller GEDs undergoing the same pressure driving as here may react more violently and collapse to smaller minimum radii leading to the absence of both the jet impact and partial deposition regimes. Therefore, a more general approach for regime identification is needed.
Although the main driver of the overall release process is the shock wave, the speed at which the drop is impacted by the surrounding liquid, 
$u_{i}$, is the actual variable dictating the drop release regime. More precisely, it is the competition between inertial, capillary and viscous forces that dictates the regime taking place. The drop impact and jet impact regimes are clearly inertia-driven, while in the partial deposition capillary effects play a much more important role. In drop impact theory, both the Weber and Reynolds numbers are leveraged to describe the competition between all three forces. Past studies have investigated the transition between the splashing and deposition of drops on flat surfaces and found that this happens for a given experimental parameter 
$K$ (Mundo, Sommerfeld & Tropea Reference Mundo, Sommerfeld and Tropea1995; Josserand & Zaleski Reference Josserand and Zaleski2003; Zhang, Zhang & Zheng Reference Zhang, Zhang and Zheng2008; Cimpeanu & Moore Reference Cimpeanu and Moore2018; Sykes et al. Reference Sykes, Fudge, Quetzeri-Santiago, Castrejón-Pita and Castrejón-Pita2022). Thoroddsen et al. (Reference Thoroddsen, Thoraval, Takehara and Etoh2011) defined the splashing parameter as a local Weber number, 
$We_{l} = K^{2}$, characterising the sheet jet, which is done by expressing the ejecta sheet thickness as 
$\delta _{s} \sim \sqrt {\mu 2r_{d,0} /(\rho u_{i})}$ and the ejecta velocity as 
$u_{s} \sim u_{i}\sqrt {Re}$. By substitution, the splashing parameter becomes
\begin{equation} K = \sqrt{We_{l}} = \sqrt{\frac{\rho \delta_{s} u_{s}^{2}}{\gamma}} = We^{1/2}Re^{1/4}. \end{equation}
Here, 
$K_{c}$ is a critical number, the value of which has been found to vary in experiments between a few tens and a few hundreds (Josserand & Zaleski Reference Josserand and Zaleski2003; Cimpeanu & Moore Reference Cimpeanu and Moore2018; Sykes et al. Reference Sykes, Fudge, Quetzeri-Santiago, Castrejón-Pita and Castrejón-Pita2022). Over the wide range of experiments conducted herein, the transition between the splashing dynamics of the drop impact regime and coalescence witnessed within the partial deposition regime occurs for a critical number 
$56 < K_{c} < 96$, which is close to values reported previously (Mundo et al. Reference Mundo, Sommerfeld and Tropea1995).
Another general parameter playing a major role in defining the drop dynamics regime is the shape of the air bubble surface as it makes contact with the drop, may it be a flat interface or the tip of a jet. Here, the local curvature of the bubble surface, 
$\mathcal {K}_{b}$, is measured at the contact point with the drop and normalised to the drop's curvature, 
$\mathcal {K}_{d} = 1/r_{d,0}$. The local curvature is assessed from the inside of the bubble, meaning that a convex and concave interface contours yield positive and negative curvatures, respectively. Figure 10(b) displays the general regime map, where a clear demarcation between the different regimes is found. Clearly the Weber and Reynolds numbers play an important role in defining whether inertia will be driving strong dynamics such as the formation of a sheet jet or the deformation of the drop upon impact with the jet, as observed by the vertical limit set at 
$K_{c}=76$. The difference between the drop impact and jet impact regimes can be explained by the local curvature of the bubble interface-to-drop curvature ratio. A value larger than one (
$\mathcal {K}_{b}/\mathcal {K}_{d} > 1$) suggests that the local radius of jet impacting the drop is smaller than the drop radius, inducing the jet impact regime, for which the drop adopts a bag-like shape. For values smaller than zero (
$\mathcal {K}_{b}/\mathcal {K}_{d} < 0$), where zero means that the bubble interface is flat and negative values means the interface has not yet started becoming convex, the drop is released in the drop impact regime. Between zero and one (
$0 < \mathcal {K}_{b}/\mathcal {K}_{d} < 1$), the bubble interface is convex and larger than the drop. The impact of the curved interface with the drop causes the drop to deform in a flat bag-like shape Sykes et al. (Reference Sykes, Fudge, Quetzeri-Santiago, Castrejón-Pita and Castrejón-Pita2022). This limit region is however hardly distinguishable in the radiographs at current resolutions. There are no experimental points in the top left side of the regime map, where jets must impact the drop at low speeds to enable coalescence between the two, which has not been witnessed experimentally.
3.5. Drop release
Finally, the implications of the identified regimes for shock-driven GED dynamics to the release of the drop within the surrounding liquid are investigated. The drop release is investigated numerically by tracking the volume of liquid initially within the drop, which can be traced by the colour applied to it. The assumption of axisymmetry in the numerical simulation is likely to influence the problem. Nonetheless, the good agreement between the experimental radiographs and numerical simulations so far suggests that the general dynamical behaviour of the drop release within the surrounding liquid must be well captured by the numerical simulation. Since the drop dissemination is a relatively slow process, the effects of an impulse on the drop dynamics is investigated. The impulse effectively comprises the combined effects of a laser-induced shock wave, the growth and collapse of a laser-induced vapour bubble and a collapse-induced shock wave.
Figures 11(a)–11(c) display the dissemination and mixing through time of drops in the drop impact,partial deposition and jet impact regimes, respectively. All three drops have an initial radius of 0.44 mm and are within bubbles of 1.2 mm radius. The shock driving conditions are selected to induce impact speeds close or equal to those within the examples presented in all three regimes. In each case, the shock wave comes from the left and the initial location of the drop and bubble are the same as the ones used for the test case displayed in figure 3 for the drop impact and partial deposition regimes, and the impact speeds are 
$u_{i}= 10.98\ \textrm {m} \textrm {s}^{-1}$ and 
$u_{i}= 2.94\ \textrm {m}\ \textrm {s}^{-1}$, respectively. In the jet impact regime, the drop and bubble are initialised at the same location as the test case displayed in figure 7 and the jet impacts the drop at 
$u_{i}= 14.31\ \textrm {m}\ \textrm {s}^{-1}$.
Figure 11. Temporal evolution of the dissemination of a drop's volume fraction within the (a) drop impact (
$We_{d} = 1940$, 
$Re_{d} = 11\ 000$), (b) partial deposition (
$We_{d} = 104$, 
$Re_{d} = 2580$) and (c) jet impact (
$We_{d} = 2460$, 
$We_{j} = 842$, 
$Re_{d} = 12500$, 
$Re_{j} = 4280$) regimes for dimensionless times of 
$u_{i}t/r_{d,0} = 0, 4, 8, 12, 16$ and 20. The initial location of the drop is indicated as a dashed circle on every frame. (d) Time evolution of the dissemination of the drop's mean density of distribution, 
$D$, within the surrounding liquid normalised to its initial state, 
$D_{0}$. The colour bars display the maximum spreading radius of the drop contents, normalised to the drop's initial radius.
The drop dissemination of the drop impact regime displayed in figure 11(a) is quantified by measuring its mean density of distribution (Stone & Stone Reference Stone and Stone2005) normalised to its initial value before contacting the surrounding liquid, 
$D/D_{0}$. The mean density of distribution is computed at each time step by averaging the volume fraction of drop contents by grid cell over the region of interest, while only considering cells for which the colour associated with the drop is non-zero. The time evolution of 
$D/D_{0}$ and dimensionless maximum spreading radius of drop contents, 
$r_{D,max}/r_{d,0}$, are shown in figure 11(d). Once the shock wave initiates the bubble collapse and the surrounding liquid impacts the drop, at 
$u_{i}t/r_{d,0}=0$, the major part of the drop quickly coalesces at the proximal bubble side. While splashing unfolds, drop contents are entrained to the distal bubble side by the sheet jet. In the second frame of figure 11(a), the mean density of distribution of the drop has already dropped to 
$0.46D_{0}$ thanks to the strong inertial driving of the splash. The splash transports the drop contents across the bubble, rapidly increasing the maximum spreading radius to 
$r_{D, max} \approx 5r_{d,0}$. Dissemination of the drop is further enhanced by the central jet that forms at the latest stage of splashing and travels both upstream and downstream. It deforms the main bulk of the drop upon impact and brings the mean density of distribution of the drop down to a value of 
$0.35D_{0}$. At 
$u_{i}t/r_{d,0}\approx 12$, the drop is impacted by the second shock wave emitted at the collapse of the vapour bubble. This induces a second contraction of the GED which briefly reconcentrates part of the drop contents, driving the mean density of distribution of the drop up to 
$0.40D_{0}$. Once the second shock wave has passed, the bubble relaxes and the main bulk of the drop is entrained by the vapour bubble-induced flow towards the origin of the shock wave. This increases the distance over which the drop disseminates to almost 12 times the initial drop radius. Parts of the drop contents can be seen clustering to the right, at the gas–liquid interface, on the two last frames of figure 11(a). At that point, the impulsive driving is over and the mean density of distribution of the drop has reached 
$0.29D_{0}$.
The partial deposition regime, displayed in figure 11(b), shows a much milder dissemination rate in figure 11(d). The bubble collapse is driven by a weaker impulse, which results in coalescence of the drop without splashing at the proximal bubble side, as soon as the bubble surface impacts the drop. The post-collapse expansion of the bubble interface induces stretching and necking of part of the drop. In the third and fourth frames of figure 11(b), one observes the formation of a daughter droplet. However, before the daughter droplet can detach from the main bulk of the drop, the second collapse of the bubble induces the total coalescence of the drop. Meanwhile, the main bulk of the drop is entrained by the vapour bubble-induced flow towards the shock wave origin, as observed in the last two frames of figure 11(b). Quantitatively, the dissemination of the drop is clearly not very efficient as the mean density of distribution of the drop slowly diminishes to 
$0.80D_{0}$ and the maximum spreading radius of the payload barely reaches 
$r_{D, max}= 2r_{d,0}$ at the end of the impulsive driving.
In the jet impact regime displayed in figure 11(c), the release dynamics of the drop are quite different from the ones previously introduced. The impact of the jet with the drop first deforms the drop before setting the entire bulk in motion at 
$u_{i}t/r_{d,0}\approx 8$. At this instant, the maximum spreading radius of the drop increases thanks to the liquid being pushed to the sides and the drop adopting a bag-like shape. As time passes, the bag becomes larger and thinner, inducing a steady dissemination of the drop until the mean density of distribution of the drop reaches 
$0.45D_{0}$ at 
$u_{i}t/r_{d,0}\approx 12$. The bag-shaped drop contacts the outer liquid during the second air bubble collapse. The increase in the mean density of distribution of the drop due to the bubble collapse-induced drop focusing is not as appreciable as in the drop impact case. The second shock wave, however, further improves the mixing of the drop to the surrounding, bringing the mean density of distribution of the drop down to 
$0.29D_{0}$. At that point, the impulsive release is over and a large amount of drop contents can be seen clustering at gas–liquid interfaces in the last frame of the image sequence of figure 11(c). This regime suggests a better distribution of the payload, although the maximum spreading radius of the drop is limited to approximately seven initial drop radii.
One observes a major difference in the mixing efficiency between the partial deposition regime and the two other regimes. The dissemination of the drop is clearly enhanced by inertia, and, unsurprisingly, higher Weber and Reynolds numbers provide a better framework for mixing in general. In the drop impact and jet impact regimes, the provided impulse is similar and the dissemination, although slower in the jet impact regime, ends up being the same. The advantage of the jet impact regime is that the distribution of the payload appears qualitatively more uniform and closer to the initial location of the GED, while in the drop impact regime most of it is entrained upstream towards the shock origin.
3.6. Discussion
Many applications involving fluidic transport, such as targeted drug delivery, would benefit from downsizing the GED significantly to the micrometric scale. At such a scale, non-negligible capillary effects on the GED release dynamics are expected due to the significant increase in interface curvatures and Laplace pressure and the required presence of stabilising surfactants to prevent bubble dissolution. Possible effects might be: a larger value for the critical number, 
$K_{c}$, which describes the transition from the partial deposition to the drop impact regime; or an increased energy requirement for the generation of jets and sheet jets such as in the jet impact and drop impact regimes, respectively. The GED dynamics are also expected to be much faster, increasing the difficulty to observe them experimentally. Nonetheless, here, the study of millimetric GED offers a first experimental framework in which the dynamics can be clearly assessed and constitutes a first step in understanding the shock-induced release dynamics of GEDs loaded with a single drop having a drop-to-bubble size ratio of 
$0.2 \leq \varepsilon \leq 0.4$.
In practice, stabilised micrometric GEDs can be loaded with single or multiple stabilised drops having a larger drop-to-bubble size ratio range (Poortinga Reference Poortinga2013; Moreno-Gomez et al. Reference Moreno-Gomez, Athanassiadis, Poortinga and Fischer2023). Although investigating GEDs that carry single drops is relatively straightforward, the presence of additional drops complicates the problem significantly, in most cases inducing a loss of axisymmetry. This is problematic in experiments, where depth perception is absent because X-ray phase contrast imaging integrates all phase discontinuities along the beam, and in the numerics where computationally costly three-dimensional simulations would be required. Three-dimensional, complex drop–drop interactions within the bubbles following the shock wave impact on the GEDs are also expected to be arduous to model analytically.
In the present work, the effects of the first laser-induced shock wave cannot be completely decoupled from the effects of the vapour bubble growth and collapse and the second vapour bubble collapse-induced shock wave. Decoupling single shock–GED interactions would be more reflective of the physics taking place in applications involving micrometric GEDs, where the dynamics are close to the time scale of the shock wave, and therapeutic shock wave sources, which often emit single pressure pulses. Despite these limitations, the effects of shock waves on the GED are expected to be dominant with respect to the velocity and pressure field caused by the vapour bubble dynamics. Indeed, the bubble compression phase and jet observed in figure 7 are typical of shock-induced bubble collapses and their time scale (Johnsen & Colonius Reference Johnsen and Colonius2009), with the jet attaining a velocity an order of magnitude larger than the vapour bubble-induced radial flow field. Therefore, the dynamics of a GED interacting with an isolated shock wave are not expected to be fundamentally different from those observed herein. The release and mixing of the drop should be the only considerably affected metric, evolving at a slower pace without the contribution of a second shock wave.
The drop release is solely investigated numerically due to experimental limitations. Indeed, eosin, which emits light upon contact with laser light, has unsuccessfully been tested as GED drop composition. The leaking of the fluorescent substance to the surroundings already at the GED production stage made it impossible to distinguish the drop from the fluorescent material in the surroundings. The use of stabilised fluorescent loaded antibubbles or GEDs could overcome this limitation.
4. Conclusion
The impulsive release of GEDs has been investigated through high-speed X-ray phase contrast imaging and numerical simulations. The impulse in this experiment is mainly provided by the combination of a relatively thin laser-induced shock wave, the radial flow induced by the expansion of a vapour bubble and a second shock wave emitted at the vapour bubble collapse. The undistorted optical access to all phase discontinuities along the X-ray path unravels the drop dynamics within the collapsing bubble. The radiographs are leveraged to corroborate the numerical method, with an excellent agreement found between the two.
Three different drop release regimes are identified for the first time, namely the drop impact, partial deposition and jet impact regimes. The drop impact and partial deposition regimes display strong similarities with the canonical configuration of a drop impacting a deep liquid pool at high and low Weber and Reynolds numbers, respectively. Existing scaling laws are successfully applied to describe the spreading of the necking region between the drop and the bubble surface in the first instants of contact. In the drop impact regime, a sheet jet forms and evolves into a crown within the bubble. This phenomenon promotes the dispersion of a small part of the drop and, overall, high Weber and Reynolds numbers promote the dissemination and mixing of the drop.
In the partial deposition regime, only parts of the drop coalesce with the surrounding liquid after the passage of the first shock wave. Strong similarities are reported with the partial coalescence of drops encountered in low-velocity drop impact experiments. However, here, the formation of a daughter droplet is driven by the oscillating motion of the bubble interface and not by capillary waves. In most experiments, the drop stretches but does not detach to form a daughter droplet within the bubble because of the second contraction phase of the bubble. Although the second shock wave enhances the mixing of the drop with the surrounding liquid, at that stage, the mixing is less effective by at least 76 % in comparison to an inertia-driven processes such as those in the drop impact regime.
In the jet impact regime, the jet splashes onto the liquid drop, showing similarities with binary drop collision dynamics. In the present work, a fitted scaling law, depending on both the drop characteristic time and drop-to-jet size ratio, is proposed for the displacement of the drop in the direction of propagation of the jet. This sheds light on the dynamics involved in the injection of the drop into the surrounding liquid. Although the mixing rate is slower than in the drop impact regime, at latter times the mean density of distribution of the drop is equivalent. The different regimes are found to be solely dependent on the Weber and Reynolds numbers and shape of the bubble interface at the initial point of contact with the drop. For a critical number 
$56 < K_{c} < 96$, the partial deposition regime transitions to inertia-driven regimes. The value found for the critical number is in the same order of magnitude as values reported previously for drops impacting flat liquid surfaces. If the local curvature of the bubble interface is higher than the drop curvature, a jet, smaller than the drop, impacts the drop leading to dynamics within the jet impact regime. For values lower than the drop curvature, the drop effectively impacts a larger surface and splashing or coalescing dynamics take place.
The present work lays the ground for understanding and controlling instantaneous targeted-payload release from GEDs using impulsive stimuli such as sound waves.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.1124.
Funding
The authors acknowledge the financial support from ETH Zurich, the ETH Zurich Postdoctoral Fellowship program and access to beamtime at beamline ID19 of the ESRF in the framework of the Shock BAG MI1397. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.
Declaration of interest
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available under reasonable request.
Appendix A
A grid convergence study is performed on the same GED as displayed in figure 3. However, since the initial conditions inducing the shock wave and growth of the vapour bubble are grid-dependent and would require an unrealistically fine grid to converge, the convergence of the GED studied herein is investigated using a plane, sustained shock wave of post-shock pressure 
$p_{2} = 20$ MPa, stronger than an impulsive shock wave. The shock driving induces a faster sheet jet which has formed at 
$t/(r_{b,0}\sqrt {\rho /p_{b,0}})=0.05$, much faster than in figure 3. Figures 12(a) and 12(b) show the volume fraction of water, 
$\alpha _{w}$, and speed along the 
$x$-axis, 
$u$, at the location, 
$x$, of the GED within the region of interest, respectively. Quantities are extracted along the horizontal dotted black line located at 
$y=r_{0}$ and indicated in the inset of figure 12(b). Figure 12(c) shows the contour line for 
$\alpha _{w} = 0.5$. Small discrepancies are visible at the sheet jet tip, which is much thinner than the droplet radius. However, these discrepancies are deemed well within the experimental uncertainty. The solution appears to have reasonably converged for 104 grid points per drop diameter. A small discrepancy can be observed for 
$x=301.63$ mm where the dotted line crosses the bubble cusp. In addition, some fluctuations in the velocity field can be observed at that location and within the gas phase of the GED shown in blue in the inset of figure 12(b). At these locations, a finer grid could improve the solution. However, this would be much more computationally demanding without any noticeable benefits to the present work. Therefore, given that a good agreement between experiments and numerics is already found, and that numerous simulations need to be carried out, the resolution of 
$n=104$ is deemed appropriate for the present study.
Figure 12. Grid convergence study indicating (a) the volume fraction of water, 
$\alpha _{w}$, and (b) the flow speed, 
$u$, along the 
$x$-axis. Theses quantities are probed one drop diameter away from the GED axis, indicated by the dark dotted line in the inset of (b). (c) Contour lines of the GED for 
$\alpha _{w} = 0.5$. Results are displayed for different numbers of grid points per drop diameter, 
$n$. A good compromise between accuracy and computational load is found for 
$n=104$.
Appendix B
The initial conditions are validated in figure 13 by comparing the size evolution of the vapour bubble and gas phase of the GEDs over time between the simulation and corresponding experimental case. The metric is the equivalent radius for all and is extracted from their area. The relative deviation, 
$\epsilon$, is mostly kept within 7 % for 
$t/\tau _{c} < 0.38$, with a 11 % spike for the vapour bubble because of its relatively small size at early times. Beyond 
$t/\tau _{c} = 0.4$, the relative deviation on the vapour bubble's radius starts increasing as the simulation diverges from the experiment as discussed in § 2.2.3. The divergence between the simulated and experimentally recorded GED dynamics is also visible and appears at an earlier time. The over-expansion of the air bubble for the simulations can causes the relative deviation to grow to 10 % for the GED observed in figure 7, and almost 20 % for the GED observed in figure 3. Although this is a considerable divergence, considering the imperfect experimental initial conditions with respect to the simulation, the evolution of the deviation through time is considered herein acceptable.
Figure 13. Relative deviation of the simulated vapour bubble and GEDs equivalent radius with respect to the experimental recordings. The shaded area shows the minimal experimental uncertainty on the radius of one pixel.
Open access
