Rayleigh–Taylor instability of cylindrical water droplet induced by laser-produced cavitation bubble

Abstract This study gives insights into the interfacial instability of water droplets through a combination of laboratory experiments, numerical simulations and analytical modelling. An experiment is conducted in a narrow gap between two plates to model the two-dimensional cylindrical geometry, and a pulsed laser beam is focused inside a water droplet to generate a cavitation bubble. Three distinct characteristics of droplet deformation can be observed: (i) splashing; (ii) ventilating; and (iii) a stable state. In addition, an analytical model considering the Rayleigh–Taylor instability and bubble oscillation is developed based on the assumption that fluid is inviscid and incompressible. The analytical model is solved to obtain a phase diagram describing three distinct phenomena. Two dimensionless parameters, $\hat {\eta }_{ins}$ and $\hat {\eta }_{sta}$, are used to determine the boundaries between different regimes. The parameter $\hat {\eta }_{ins}$ is defined by the ratio of the perturbation amplitude to the difference between the droplet and bubble radii, whereas the other parameter $\hat {\eta }_{sta}$ is defined by the ratio of the perturbation amplitude to the initial droplet radius. Interfacial instability is induced when the perturbed droplet surface penetrates the bubble boundary, namely, $\rvert \hat {\eta }_{ins} \rvert \geq 1$. Splashing and ventilating phenomena occur for $\hat {\eta }_{ins}\geq 1$ and $\hat {\eta }_{ins}\leq -1$, respectively. A stable state occurs when droplet fragmentation does not appear despite small-amplitude perturbations for $\hat {\eta }_{sta}\leq 0.1$. There is a transition zone between the ventilating and stable state, which is bounded by $\hat {\eta }_{ins} = -1$ and $\hat {\eta }_{sta} = 0.1$. Finally, the phase diagram is verified by the experimental and numerical results.


Introduction
The Rayleigh-Taylor (RT) instability appears in various natural and industrial processes, such as inertial confinement fusion (Boudesocque-Dubois, Gauthier & Clarisse 2008; Kuhl et al. 2013), supernova collapse (Plewa 2007;Hicks 2015), and ocean dynamics (Veron et al. 2012;Veron 2015). When a heavier fluid lies above a lighter fluid, the system is in a gravitationally unstable state, and any small perturbation of the interface will grow with time, leading to the RT instability. The spherical RT instability is also referred to as the Bell-Plesset effect, which accounts for the influence of global convergence or divergence of fluid on the growth of surface perturbations. Bell (1951) derived the governing equations to study the perturbation growth induced by the radial motion of the interface in spherical and cylindrical geometries assuming inviscid, incompressible and immiscible fluid flowing irrotationally. Later, Plesset (1954) improved the model and obtained differential equations for the perturbation when the amplitude was small compared to its wavelength. Thereafter, Bell's results were extended to an arbitrary number of fluid layers in planar, cylindrical and spherical geometries, while considering the effects of liquid compressibility, viscosity and surface tension (Prosperetti 1977;Mikaelian 2005;Yu & Livescu 2008;Forbes 2011;Bakhsh et al. 2016;Zhou 2017).
The spherical RT instability of a water droplet induced by internal volume oscillation has attracted much attention and is very important for industrial applications such as premixed combustion (Chertkov, Lebedev & Vladimirova 2009), underwater explosion (Geers & Hunter 2002) and industrial coating with thin liquid films (Livescu & Schwartz 2011). The volume oscillation inside a droplet is induced by expansion and contraction of a cavitation bubble, which is commonly generated by a pulsed laser and electric discharge. For the absorbed energy of a droplet, the transition pathways include plasma formation, phase changes and subsequent mechanical behaviours such as shock-wave emission, bubble oscillation and droplet deformation. The characteristics of the generated plasma have been investigated regarding the propagation velocity, plume shape and density (Eickmans, Hsieh & Chang 1987;Hsieh et al. 1987;Reijers, Snoeijer & Gelderblom 2017). The subsequent expansion of plasma generates shock waves propagating around the droplet and a cavitation bubble oscillating inside the droplet (Kafalas & Ferdinand 1973;Kafalas & Herrmann 1973;Stan et al. 2016a,b). Obreschkow et al. (2006) investigated the dynamics of a spark-generated cavitation bubble in a water droplet. By generating nearly spherical droplets under microgravity (Kobel et al. 2009), they found that (i) two liquid jets are generated and then escape from the droplet when the cavitation bubble collapses toroidally, (ii) the bubble's lifetime in the droplet is shorter than that in an unlimited environment and (iii) secondary cavitation occurs near the droplet surface owing to the reflection of the shock wave. Thoroddsen et al. (2009) and Heijnen et al. (2009) focused a laser beam near the surface of a hemispherical droplet to investigate the disruption of its surface and obtained the evolution of spray and liquid jets under different conditions. Gonzalez Avila & Ohl (2016) analysed the dynamics and fragmentation of a suspended droplet caused by a laser-induced cavitation bubble. The experimental results revealed three distinct fragmentation regimes, namely, atomisation, and (c) a stable state. The droplet and bubble surfaces are indicated using white and red arrows, respectively. A ventilating phenomenon is clearly observed in the experiment for the first time, in which the three-dimensional configuration of a spherical droplet prevents obvious interior observations. The splashing phenomenon and stable state of the droplet have also been observed for a spherical droplet.
sheet formation and coarse fragmentation. The dynamics of a free-falling drop impact by a pulsed laser were investigated, and the experimental results identified two RT instabilities of different origins as the prime cause of fragmentation (Klein et al. 2015(Klein et al. , 2020. The aforementioned studies show that the deformation of water droplets is usually characterised by vaporisation, liquid jets and water films. Recently, Zeng et al. (2018) experimentally and numerically found that liquid jets are excited by spherical-RT instability. The jet is generated when a pressure impulse focuses on the troughs of the perturbed surface induced by the interfacial instability. A critical condition determining the RT instability is predicted using an indirect method, such as the observation of the generation of jets, since the perturbation amplitude at the corrugated surface is hardly measured in experiments for a spherical droplet. Quantitative measurements of surface perturbations are helpful for understanding the dynamics of water droplets.
In the present study, we attempt to use a two-dimensional (2-D) cylindrical geometry to provide clearer experimental observations for the interfacial instability of water droplets. Experiments are conducted in a narrow gap between two plates, and three distinct characteristics of droplet deformation are obtained: (i) splashing; (ii) ventilating; and (iii) a stable state (see figure 1). In figure 1(a), the droplet fragments when the bubble boundary touches the disturbed droplet surface and then the ejections splash quickly owing to the inertia of the bubble expansion. When the perturbed droplet surface penetrates the bubble boundary during bubble contraction, the ambient air enters the bubble, resulting in a ventilating phenomenon owing to the pressure difference. Large-amplitude perturbations are observed in the cases of splashing and ventilating. The stable state occurs when the droplet surface maintains stable despite small-amplitude perturbations. Consequently, these perturbations in different phenomena are clearly observed by using a 2-D cylindrical geometry, so that they can be accurately measured and quantified in the present study. We aim to distinguish the regimes of stability and instability for a droplet. The structure of the present study is shown as follows: § 2 describes the methodologies of the experiment, numerical simulation and analytical modelling. The corresponding results and discussion of the three approaches are given in § § 3.1, 3.2 and 3.3, respectively. The phase diagram of interfacial instability is evaluated in § 3.4. The work is summarised in § 4.  Figure 2. Experimental set-up for observing the dynamic behaviours of 2-D cylindrical water droplets caused by laser-induced cavitation bubbles: (a) top view and (b) side view. The red and yellow light paths are from the pulsed laser and LED lamp, respectively. Mirrors 1 and 2 are used to reflect the laser light and the LED light, respectively.

Experiments
In the experiments, droplets of various sizes are generated by injecting water in a narrow gap between two 10 mm thick acrylic plates. The size of the narrow gap is 200 mm long, 200 mm wide and 0.8 mm deep and placed horizontally on a test platform, as shown in figure 2. To generate a cavitation bubble, a pulsed laser (Nd:YAG, λ = 1064 nm, 10 ns pulse width; Spectra-Physics Ltd., USA) passes through the lower wall of the gap and focuses on a square piece of aluminium film attached to the wall. This is done through a convex lens with a focal distance of 665 mm and Mirror 1 with a diameter of 25 mm, as shown in figure 2(a). The maximum energy of the laser is approximately 2.0 J, as measured by an energy meter, and the diameter of the focal spot of the laser is approximately 2 mm. To observe the top view of the droplet, two mirrors (namely Mirror 2 with dimensions of 150 mm × 100 mm) are used to change the light path so that a high-speed camera (Phantom v2512; Vision Research Ltd., USA) and light-emitting diode (LED) pulsed lamp (Constellation 120, Veritas, USA) are placed horizontally. The high-speed camera coordinates with a synchroniser (Berkeley Neceonics Corporation, BNC575-8) when the laser is triggered. The laser path is inclined slightly toward the LED lamp to prevent the laser beam from entering the high-speed camera. As shown in figure 2(b), the LED lamp and high-speed camera are placed at the same height as the narrow gap while the side view is observed.

Numerical simulations
The OpenFOAM framework is used to perform numerical simulations to analyse the interfacial instability of a water droplet. The volume of fluid (VOF) and large-eddy simulation (LES) methods are applied to capture the air-water interface and small-scale flow structures (Suponitsky et al. 2017;Wang et al. 2017;Zeng et al. 2018). In the simulations, the gas and liquid phases are treated as compressible Newtonian fluids while considering heat transfer. The governing equations for continuity, momentum and energy are where ∇ is the gradient operator, t is the time, U is the velocity vector, ρ is the mixture density, ρ = ρ liquid α + ρ gas (1 − α), α is the volume fraction such that α = 1 for the liquid and α = 0 for the gas, p is the pressure, μ is the mixture dynamic viscosity, μ = αμ liquid + (1 − α)μ gas , μ liquid and μ gas are the dynamic viscosities for the liquid and gas phases, respectively, g is the acceleration of gravity, σ is the surface tension, κ is the curvature of the interface, T is the temperature, C v is the mixture specific heat capacity, C v = C v.liquid α + C v.gas (1 − α), C v.liquid and C v.gas are the specific heat capacities for the liquid and gas phases, respectively, λ is the mixture heat conductivity, λ = λ liquid α + λ gas (1 − α), λ liquid and λ gas are the specific heat capacities for the liquid and gas phases, respectively and k is the kinetic energy per unit mass, k =|U| 2 /2. The pressure-based implicit splitting of operators (PISO) algorithm embedded in the multiphase solver compressibleInterFoam is used to solve this transient flow problem assuming that the two phases are homogeneous (Issa 1986). The transport equation of the phase fraction is where the relative velocity U r is defined as U r = U liquid − U gas . The third term in (2.4) is an artificial compression term that sharpens the interface between the two phases, and Ur is defined as where c α is the artificial compression coefficient. Finally, the equation of state (EOS) for either phase is given by where ρ 0.phase and p 0.phase are the initial density and pressure, respectively, of the given phases and ψ is the compressibility coefficient.
Favre filtering is used to accurately predict large-scale turbulent eddies. Hence, the transport equations of continuity, momentum, energy and volume-fraction are filtered as where the tilde-annotated quantities are the Favre (density weighted) filtered quantities, the overbar quantities are the LES physical space (unweighted) filtered quantities, the subscript SGS represents the subgrid scale, is the stress tensor, H is the heat flux and β is the mass flux. The one-equation eddy-viscosity model is used for the subgrid-scale modelling, and all the partial differential equations are solved using the finite-volume method. In summary, the governing equations are solved through the following steps (Ye et al. 2019): (i) solve the transport equation in (2.10) to obtain α; (ii) solve the filtered continuity equation in (2.7), the momentum equation in (2.8) and the energy equation in (2.9) in sequence; (iii) fix the pressure and update ρ phase through the EOS in (2.6); (iv) execute the pressure-correction loop until the accuracy requirement is satisfied; (v) fix the temperature and update ρ phase through the EOS in (2.6); (vi) update the mixture density ρ; (vii) end the pressure-correction loop and PISO loop, then advance the time step.
The implicit Euler scheme is used to discretise the temporal derivatives. A second-order van Leer scheme and a second-order upwind scheme are used to discretise the convective terms in the volume-fraction equation and the momentum equation, respectively. In the simulations, the other coefficients are as follows: the density of water ρ 0.liquid is 998 kg m −3 ; dynamic viscosity μ is 1.307 × 10 −3 Pa s; gravity acceleration g is 9.81 m s −2 ; surface tension σ is 0.0728 N m −1 ; specific heat capacities C v are 4199 J (kg K) −1 for water and 723 J (kg K) −1 for gas; and heat conductivities λ are 0.599 W (m K) −1 for water and 0.034 W (m K) −1 for gas.
The computational domain in the numerical simulations is shown in figure 3. A quarter of the water droplet is simulated to improve the calculation precision and computational efficiency. The size of the computational domain is 15 mm (L) by 15 mm (W) by 0.8 mm (D). The initial cavitation bubble is regarded to be a gas sphere with high internal pressure at the centre of the narrow gap in the simulations. The O-block structured grid is used, and the total number of cells is approximately 916 800. The detailed grid structure and mesh are shown in figure 3(b) and 3(c), respectively. The cell sizes are initially set to approximately 0.01 and 0.03 mm in the radial direction near the bubble and droplet, respectively. To precisely capture the air-water interface, the grids are denser near the walls as well as in  Figure 4 shows a schematic of the interfacial evolution of a water droplet driven by the oscillation of a cavitation bubble. The droplet surface is initially perturbed by the reflection of shock waves. As the cavitation bubble oscillates, the surface perturbation grows rapidly when the droplet surface is accelerated by ambient air. A multi-mode sine function is used to describe the surface perturbation, and we denote by η the amplitude of the crest. Interfacial instability occurs when the perturbed surface touches the bubble boundary and the droplet fragments. In other words, the amplitude of the perturbation is equal to the difference between the bubble radius and droplet radius.

Analytical model for interfacial instability
To investigate the evolution of the perturbation at the droplet surface, we develop an analytical model based on the assumption of an inviscid and incompressible fluid in a 2-D cylindrical geometry. It comprises (i) a small-amplitude model to describe the perturbation growth and (ii) a bubble-oscillation model to describe the dynamic characteristics of the droplet. In the small-amplitude model, the Bell equation (Bell 1951) is used for a 2-D cylindrical geometry, such that where η(t) and R d (t) are the perturbation amplitude and the droplet radius, respectively (the dots over these terms denote differentiation with respect to time), n is the mode number, that is, the number of wave crests at the droplet surface (see figure 4) and A is the Atwood number. The Atwood number is defined as where ρ 1 and ρ 2 are the densities of the ambient air and the droplet, respectively. An equation governing the evolution of the droplet radius is needed to close the system.  Figure 4. Schematic of the interfacial evolution of water droplets driven by the oscillation of cavitation bubbles. In the figure, η is the amplitude of surface perturbation, R d is the droplet radius, R b is the bubble radius, P R is the pressure on the bubble boundary, P d is the pressure on the droplet surface and n is the mode number, i.e. the number of wave crests at the droplet surface.
Ignoring the compressibility of the liquid, we can obtain the size and velocity of the droplet surface assuming a constant droplet volume, namely, where R is the radius and the subscripts b and d correspond to the bubble and the droplet, respectively. Next, a model for bubble oscillation is developed. We express the governing equations in cylindrical geometry as where p is the liquid pressure, ρ is the liquid density, u is the liquid velocity in the radial direction and t is the time. The compressibility of the liquid is neglected and (2.16) is integrated with respect to time. Therefore, the oscillation of a bubble in a water droplet is modelled asṘ where P R is the pressure on the bubble boundary, P d is the pressure on the droplet surface, R d is the droplet radius and R b is the bubble radius. We express P d as where P air is the pressure of the ambient air (1.01 × 10 5 Pa herein), σ is the bubble surface tension coefficient and μ is the dynamic viscosity. As mentioned previously, a narrow gap is used in the experiments. Consequently, viscous boundary layers arise between the walls, leading to a viscous force. The viscous force is estimated and added to the analytical model. Hence, the pressure on the droplet surface is described as is the surface area of the droplet in the radial direction and A 2 = 2πR d h is the surface area of the droplet in the direction of the gap depth h. The resistance coefficient C f is obtained by the Blasius equation: (2.20) In addition, we express P R as where P in is the pressure inside the bubble. For an ideal gas, the EOS inside the bubble is where v and T in are the molar volume and internal temperature of the bubble, respectively where N A = 6.02 × 10 23 is the Avogadro's number, V is the bubble volume and N Tot is the total number of gas molecules inside the bubble. Considering the effect of thermal conduction at the bubble boundary, the energy balance for the gas inside the bubble is given by (2.24) according to the first law of thermodynamics, where the bubble radius varies from R b to R b + ΔR b and the gas temperature changes by ΔT during time Δt: ( 2.24) where C vg is the specific heat of the gas at constant volume and S b is the area of the bubble surface. The term on the left-hand side is the variation in the internal energy of the gas during time Δt, the second term on the right-hand side is the work done by the pressure on the bubble surface and ΔQ is the heat released from the bubble to the liquid through the thermal boundary layer, which exists inside the bubble because the density and specific heat of water are much larger than those of the gas (Toegel et al. 2000;Toegel & Lohse 2003;Yasui & Kato 2012).
In the model, we assume that the temperature is spatially uniform within the bubble but varies linearly within the boundary layer. Therefore, according to Fourier's law, ΔQ can be written as where λ g is the thermal conductivity, δ is the thickness of the thermal boundary layer and T f is the temperature at infinity. The thickness δ can be obtained from the Rayleigh-Plesset Weber number P d0 R d0 σ 10 3 -10 4 ≥1.4 × 10 3 ≥1.4 × 10 3 Reynolds number R d0 √ P d0 ρ d0 μ 10 4 -10 5 ≥0.7 × 10 4 ≥0.7 × 10 4 Mode number n --≥1 is the thermal diffusivity, C p is the heat capacity at constant pressure, R bmax is the maximum bubble radius, ΔP is the pressure difference between the inside and outside of the bubble when it reaches R bmax and ρ g is the gas density, which varies with the temperature and pressure inside the bubble. In the model, the effects of evaporation and condensation inside the bubble are ignored because the rates of the phase change are slower than the surface velocity of the bubble surface (Brian & Andrew 2000). Finally, (2.11) and (2.17) form a complete system that is solved using the fourth-order Runge-Kutta-Gill method.

Analysis of the controlling dimensionless parameters
As mentioned previously, the amplitude of the perturbation is used to determine whether interfacial instability occurs at the droplet surface. It is important to discuss how controlling parameters affect the perturbation growth under different initial conditions of the cavitation bubble and water droplet. Table 1 lists the corresponding dimensionless parameters, the subscripts b and d represent the bubble and droplet, respectively and the subscript 0 represents the initial status.
First, the initial conditions of the cavitation bubble involve the initial pressure, radius and temperature. In the experiment, the initial conditions of the bubble are adjusted by changing the area of the aluminium film with constant laser energy. In this way, the initial temperature of thermal products inside the cavitation bubble is maintained since the energy of the laser remains in the experiments. The initial size of the cavitation bubble is of the order of micrometres. It is reasonable to maintain a constant initial radius and adjust the internal pressure by changing the area of the film. It is difficult to measure directly the initial status of the cavitation bubble and droplet. In the present study, the initial conditions are estimated by substituting experimental results into the analytical model (see § 3.3). Consequently, the radius R b0 is 350 μm and the pressures are 10 and 7 MP, respectively, when we choose aluminium films larger and smaller than the focal spot of the laser. The initial conditions of the droplet include the radius, pressure, density and initial perturbation. According to the experiments, the radius R d0 ranges from 1 to 10 mm, the pressure is P d0 = 1.01 × 10 5 Pa and the density is ρ d0 = 998 kg m −3 . As shown in figure 4, the droplet surface is initially perturbed by the reflection of the shock wave. First, the Tait equation is used to obtain the liquid density behind the shock wave and then the conservations of mass and momentum are employed to calculate the corresponding velocity. Finally, the initial perturbation is obtained as the product of the velocity and time interval. The time starts from the initial generation of the reflected wave until it reaches the droplet centre. The time period is estimated to be of the order of microseconds, which shows a good agreement with Period 1 mainly caused by the propagation of the shock wave shown in Appendix B. Referring to the curve of shock pressure attenuation measured by Wang et al. (2018b), the pressure of the shock wave is estimated to vary from 2 to 10 MPa at the droplet surface of various sizes. The corresponding velocities behind the shock wave are from 0.070 to 0.636 m s −1 . As a result, the amplitude of the initial perturbation ranges from 5.09 × 10 −7 to 10 −6 m, and the corresponding dimensionless value ranges from 10 −5 to 10 −4 .
The other coefficients used in table 2 include the density of ambient air ρ air = 1.20 kg m −3 , the dynamic viscosity μ = 1.307 × 10 −3 Pa s and the surface tension σ = 0.0728 N m −1 . The effect of surface tension can be ignored because the Weber number in table 1 is much larger than 1. The viscosity near the walls of the plates is considered in (2.20) and merely affects the dynamics of cavitation bubbles and droplets owing to a large Reynolds number. Consequently, the key parameters affecting the growth of the surface perturbation in the present analysis are the pressure ratio, radius ratio, initial perturbation and mode number. Figure 5 shows the splashing phenomenon of a 4.6 mm diameter droplet. The high-speed camera captures the images at a frame rate of 100 kfps with an exposure time of 300 ns. The black square at the centre of the water droplet is an aluminium film of size 2.0 mm × 2.0 mm, which is larger than the focal spot of the laser to absorb all the laser energy. Figures 5(1), 6(1) and 7(1) are captured just after the lasers are triggered. In figure 5(1), the strong flash indicates that the laser focuses on the aluminium film and that the cavitation bubble has been generated. The water droplet grows with the expansion of the bubble, and a ring-shaped jet (first jet) escapes from the main droplet (see figure 5 (5)). In addition, a shadow between the surfaces of the droplet and bubble is clearly observed, which indicates that the droplet surface is perturbed. Subsequently, the droplet fragments when the bubble boundary contacts the droplet surface and splashes quickly owing to the inertia of the bubble expansion. The cavitation bubble also breaks apart when contacting ambient air, and these ejections rapidly reach the escaped jet, forming a splashing phenomenon. The interested reader is referred to Appendix B for experimental and numerical investigations of the generation mechanism of the first jet shown in figure 5(5). Figure 6 shows the ventilating phenomenon of a 12 mm diameter water droplet. The settings of the high-speed camera are the same as those in figure 5. The size of the aluminium film is 3.5 mm × 3.5 mm. In figure 6(1), a bubble cloud with an annular shape is captured around the film, which can be generated by the expansion wave after the shock wave reflects on the droplet surface (Wang et al. 2018a). In figures 6(2) and 6(3), the water droplet grows with the expansion of the laser-induced bubble. Just before the bubble reaches its maximum size at approximately 330 μs, a ring-shaped jet is also observed in this case, and Appendix B shows its generation process. As the cavitation shrinks, perturbation at the droplet surface grows remarkably and penetrates the bubble boundary, forming channels between the bubble interior and ambient air. The ambient air flows toward the interior of the bubble owing to the pressure difference, resulting in a ventilating phenomenon, as shown in figure 6(7). This indicates that interfacial instability occurs during the contraction of the bubble. At approximately 800 μs, the cavitation bubble collapses, generating a pressure impulse that reverses the moving direction of the droplet surface. Figure 7 presents a stable state for a droplet with a diameter of 20 mm. The settings of the high-speed camera are the same as those in figure 5. The size of the aluminium film is 2.0 mm × 2.0 mm. In figure 7(1), a bubble cloud around the film is also observed just after the laser is triggered. As the cavitation bubble expands, the droplet grows to its maximum size at approximately 440 μs. The ring-shaped jet is observed at the initial contraction stage, and Appendix B shows its generation mechanism. Most of the droplet surface remains smooth except for the large-amplitude distortion marked by the red circle in figure 7(4). This distortion is probably owing to the imperfect initial surface of water droplet and grows with the expansion of the bubble. In figure 7(5), some corrugations are observed at the droplet surface, and they grow until the bubble rebounds at approximately 1000 μs. Interfacial instability is not induced since the perturbed surface of the droplet does not touch the bubble boundary and fragmentation is not induced. After the second collapse at approximately 1200 μs, the cavitation bubble breaks up into small bubbles. As these small bubbles expand and contract, the droplet surface is slightly perturbed. Figure 8 shows comparisons of the ventilating phenomenon obtained experimentally and numerically. When similar maximum sizes and collapse times of the bubbles in the simulation match with the experimental observation, the initial bubble radius, internal pressure and internal temperature are set to 350 μm, 16.5 MPa and 1500 K, respectively, referring to the numerical simulations in the study of Zeng et al. (2018). From the figure, it is found that the numerical results show good agreement with the experimental observation regarding the dynamics of the bubble and droplet. On this basis, a detailed numerical analysis of the evolution of the droplet surface is performed. We numerically investigate the evolution of the baroclinicity and vorticity at the droplet surface. In figure 9(a-1), a surface between the bubble and droplet surfaces is induced by the expansion wave according to its propagation velocity. This indicates that the bubble clouds with annular shapes in figures 6 and 7 are induced by the tensile area behind the expansion wave. During the later stage of bubble expansion (at 300 and 400 μs), baroclinicity is clearly induced and the droplet surface is perturbed, thereby forming some corrugations owing to the generated vorticity. It is because of this that the density gradient and pressure gradient are not aligned after the droplet surface is slightly perturbed. As the bubble shrinks, the surface perturbation grows rapidly. In figure 9(8) the perturbed droplet surface penetrates the bubble boundary. Upon the collapse of the bubble, a pressure impulse inversely transforms the distribution of the baroclinicity and vorticity, leading to the ejections in figure 6(10-12).

Analytical solution of the instability model
The initial conditions of the cavitation bubble play an important role in affecting the growth of the perturbation. It is difficult, however, to directly measure or observe in the experiments. In this work, the initial conditions of the bubble are estimated while matching the experimental results by the analytical model. Figure 10(a) shows comparisons of the diameters of the cavitation bubble and the droplet obtained experimentally and analytically (the experimental results are based on figure 6). By counting the number of surface perturbations, the mode number of n = 42 is obtained and used in the analytical model. The diamonds and squares with error bars are the diameters of the droplet and cavitation bubble, respectively, which are averaged from 20 positions of the respective surface. Referring to the initial conditions in the simulation, we can estimate the initial condition of the bubble in this case when maximum sizes of bubbles and droplets in the analytical model match with the experimental observation. Here, initial pressure inside the bubble is P in0 = 10 MPa, the radius is R b0 = 350 μm and the temperature is T b0 = 1500 K. It then follows that the pressure ratio of P in0 /P d0 = 99.0 and radius ratio of R d0 /R b0 = 17.1 for the experimental observation of figure 6. To investigate the error range of the initial conditions, we substitute the maximum and minimum of error bars into the analytical model. It is found that the maximum and minimum values of R b0 are 366 and 336 μm while that of P in0 are 10.15 and 9.35 MPa, respectively. For the other settings in the analytical model, the thermal conductivity is λ g = 0.599 W (m K) −1 , the heat capacity is C v = 4190 J (kg K) −1 , the droplet radius is R d0 = 6 mm, the density is ρ ∞ = 998 kg m −3 , the dynamic viscosity is μ = 1.307 × 10 −3 Pa s, the surface tension is σ = 72.8 × 10 −3 N m −1 , the atmospheric pressure is P 0 = 1.01 × 10 5 Pa and η 0 = 2 × 10 −7 m.
As shown in figure 10(a), the bubble shrinks faster in the experiment than that in the analytical model during the contraction of the cavitation bubble. This is because the first jet takes away part of the mass and momentum of the droplet. With the exception of the errors induced by the loss of the mass and momentum, the analytical model describes the dynamics of the cavitation bubble and the droplet well. According to the direction of the surface acceleration, an unstable period of the droplet surface is qualitatively predicted as indicated by red area in figure 10(b). Here, the direction of the bubble expansion is set as positive. It is found that the surface perturbations appear to grow exponentially at the later stage of expansion and the initial stage of contraction of the bubble. Eventually, the amplitude of the dimensionless perturbations reach the first peak, η max = η max /R d0 , 0.18-fold droplet radius at the collapse of the bubble.
Next, the analytical model is solved to investigate how the mode number n and initial perturbations η 0 affect the evolution of the surface perturbation. By using the Wentzel-Kramers-Brillouin approximation, we can obtain a function of the perturbation η, namely, η(t) ∼ η 0 exp(S n (t)). Here, S n is the growth rate of the perturbation and can be expressed as where a is the acceleration, A is the Atwood number and k w = n/R d0 is the wavenumber.
We can obtain a relation between the peak of dimensionless perturbation η max and mode number n, namely, (3.2) Figure 11(a) shows the estimation of the peak of dimensionless perturbation for a 12 mm diameter droplet under different mode numbers. The other initial conditions are the same as those in figure 10. The abscissa is the mode number n and the ordinate is the peak of the dimensionless perturbation η max . The peak of the perturbation increases as an exponential function of √ n. It is relatively easy for a large mode number to result in interfacial instability of the droplet. A coefficient of 10.82 is obtained by fitting the curve using (3.2), as shown in figure 11(a). Here, the dimensionless initial perturbation is η 0 = η 0 /R d0 = 3.33 × 10 −5 . As a result, (3.2) can be effectively used to estimate the peak amplitude of the perturbation under different mode numbers.
In addition, the estimations of the dimensionless perturbation peak are shown in figure 11(b) under different initial perturbations. The solid lines are the analytical solutions of the instability model when the dimensionless initial perturbation is 0.5, 2, 4, 10 and 20 times, η 0 = 3.33 × 10 −5 . The dashed lines are the estimations obtained with (3.2) when using the coefficient of 10.82. The estimated lines agree well with the analytical solutions. This suggests that the peak of the dimensionless perturbation is proportional to the initial perturbation and that (3.2) accurately describes how the mode number n and initial perturbation η 0 affect the perturbation growth.

Phase diagram of interfacial instability
The analytical model is solved with different initial conditions of the bubbles and droplets. Two dimensionless parameters of the surface perturbation are defined to determine the instability and stability regimes of the droplet surface. One parameter isη inŝ where η max is the maximum perturbation of the droplet surface and R d and R b are the radii of the droplet and bubble at this moment, respectively. The positive and negative signs of the termṘ b /|Ṙ b | represent the expansion and contraction of the bubble, respectively. This means that interfacial instability occurs when the surface perturbation touches the bubble boundary, i.e. the amplitude of the perturbation is equal to the difference between the bubble radius and droplet radius. Consequently, splashing and ventilating phenomena occur forη ins ≥ 1 andη ins ≤ −1, respectively. The other dimensionless parameter isη sta , which is used to determine the stability regime where the bubble boundary is not penetrated despite small-amplitude perturbations. Referring to the study of Zeng et al. (2018),η sta = |η max /R d0 | = 0.1 is applied as a threshold. Consequently, we obtain the phase diagram shown in figure 12(a). The mode number n is 42. The dimensionless pressure P in0 /P d0 and radius R d0 /R b0 are used as the ordinate and abscissa, respectively. In the phase diagram, the regimes of splashing, ventilating and a stable state are presented using red, blue and green, respectively. There is a transition zone between the ventilating phenomenon and stable state, which is bounded by the lines ofη ins = −1 andη sta = 0.1. Figure 12(b) shows typical examples in the respective regime. The variables E and S represent the results obtained from experiments and simulations, respectively. The splashing phenomenon occurs during bubble expansion in the high-energy and small-droplet regime. We obtain |η max /R d0 | = 0.45 for E-1 and |η max /R d0 | = 0.55 for S-1. For E-2 and S-2, the perturbations penetrate the bubble boundary during the initial contraction of the bubble, the values of |η max /R d0 | are 0.22 and 0.14, respectively. For E-3 and S-3, the perturbations grow until the cavitation bubble collapses, the values of |η max /R d0 | are 0.30 and 0.35, respectively. E-4 and S-4 are cases in the regime of the stable state, where the interfacial instability is not induced despite some corrugations at the droplet surface. The other results obtained experimentally and numerically are plotted in the phase diagram in the same way (see figure 12a). Here, the initial conditions of the cavitation bubble in the experiments are obtained by solving the analytical model, as shown in figure 10. We obtain P in0 /P d0 = 99.0 and R d0 /R b0 = 5.7 for figure 5, P in0 /P d0 = 99.0 and R d0 /R b0 = 17.1 for figure 6 and P in0 /P d0 = 99.0 and R d0 /R b0 = 28.6 for figure 7.
The phase diagram for a 2-D cylindrical droplet shows some similarities and differences compared with that of a spherical droplet (Gonzalez Avila & Ohl 2016). The splashing and stable state of the droplet are also observed for a spherical geometry. However, the ventilating phenomenon is observed in the experiment for the first time. The perturbed droplet surface in this case can touch or even penetrate the bubble boundary, forming channels between the bubble interior and ambient air. As Zeng et al. (2018) suggested, the corrugated surface of spherical droplet does not allow a closer look at the cavitation bubble in the experiment. Consequently, it is difficult to observe clearly the ventilating phenomenon in a spherical droplet. In the present study the channels are observed clearly and the perturbations at the surface are measured precisely so that the interfacial instability is accurately determined. From the perspective of energy transfer, interfacial instability occurs when the droplet obtains enough energy from the oscillation of the bubble.  Thus, the relationship between the pressure and radius ratios is obtained as where C n is the bubble-to-droplet energy conversion factor and depends mainly on the mode number. Considering the effect of the mode number in figure 11, C n can be expressed as C n = C 2 e ± √ n + C 3 .
(3.5) Figure 13 shows boundary curve lines between the stability and instability regimes of the droplet surface at n = 20, 25, 30, 40, 45 and 50. The area of the stability regime decreases as the mode number increases. This is because the wavenumber increases with the mode number, thereby increasing the growth rate of the perturbation according to (3.1). Consequently, the perturbations grow faster for a larger mode number. Overall, the estimates show good agreement with the analytical solutions. Therefore, (3.4) can be used to predict the stability and instability regimes for a water droplet.

Conclusion
In this paper, we investigate the interfacial instability of a water droplet induced by internal volume oscillation using a combination of experiments, numerical simulations and analytical modelling. In the experiments droplets of various sizes are generated by injecting water into a narrow gap, and a pulsed laser is focused on the interior of the droplet to generate a cavitation bubble. Three distinct characteristics of droplet deformation are observed, namely, (i) splashing, (ii) ventilating and (iii) a stable state.
The ventilating phenomenon is experimentally observed for the first time. As the bubble shrinks, the perturbed droplet surface moves toward its centre and penetrates the bubble boundary so that the bubble interior is connected to ambient air, resulting in the ventilation phenomenon. The phenomena of splashing and a stable state, which were reported for a spherical droplet, were also found in the present study.
To analyse the growth of perturbations at the droplet surface, an analytical model comprising (i) a linear model to describe the evolution of the interfacial instability and (ii) a bubble oscillation model is developed. A phase diagram distinguishing various phenomena of the droplet is obtained by solving the analytical model and is verified by the experimental and numerical results.
Two dimensionless parameters of surface perturbation,η ins andη sta , are defined to distinguish the different regimes of the phase diagram. Interfacial instability occurs when the perturbed surface penetrates the bubble boundary, i.e. |η ins | ≥ 1. Splashing and ventilating phenomena occur forη ins ≥ 1 andη ins ≤ −1, respectively. The droplet is considered to be stable as fragmentation does not occur despite small-amplitude perturbations, and the regime of the stable state isη sta ≤ 0.1, which is the other dimensionless parameter. There is a transition zone between the regimes of ventilating and stable states, which is bounded byη ins = −1 andη sta = 0.1.
The experimental observations and numerical results show that the location of the first jet escaping from the droplet depends on the initial shape of the water surface from the side view as mentioned in Appendix B. Therefore, it is of great interest to investigate the evolution of the first jet under different initial shapes of the water surface by changing the contact angles of the plates. Declaration of interests. The authors report no conflict of interest.

Appendix A. Grid independence
We analyse the grid independence by developing three types of grids, namely, coarse, medium and fine. Table 2 shows the numbers of cells and simulation times for each grid.
Here, the simulation time is the time taken to calculate the first pulsation of the cavitation bubble. For each grid, the resolution near the walls of the narrow gap was high enough to analyse the boundary layer. Figure 14 shows the temporal evolutions of the bubble radius and the position of the outer edge of the droplet during the first pulsation under different grids. The detailed information is given in table 3. Note that the initial conditions are the same as those in figure 8, namely, the radius is R b0 = 350 μm, the temperature is T 0 = 1500 K and the pressure is P in0 = 16.5 MPa. As shown in figure 14, all the results obtained with various grids are consistent regarding the evolution of the bubble radius. As to the outer edge of the droplet, the results of the medium grid are consistent with those of the fine grid, whereas there are some differences in the later stage of the first pulsation for the coarse grid. Consequently, considering computational cost and efficiency, we used the medium grid in all simulations.   Table 3. Details of the bubble radius and outer edge of the droplet for coarse, medium and fine grids.

Appendix B. Generation of the ring-shaped jet
To investigate the generation mechanism of the ring-shaped jet, we observed the side view of the narrow gap. The high-speed camera captures the images at a frame rate of 100 kfps with an exposure time of 300 ns. The resolution of the images is 1280 pixels × 192 pixels. In figure 15(1), the surface of the droplet in the narrow gap remains curved because of surface tension. As the bubble expands, the curved surface grows gradually, forms a water dome at 200 μs and transforms into a jet, as shown in figure 15(5). Compared with the top view in figure 5, we can conclude that the jet created at the middle of the droplet surface is the ring-shaped first jet. The position of the jet is slightly above the centre of the gap, which may be owing to the three-dimensional effect of initial expansion of the cavitation bubble. As the bubble rebounds, two further jets form at both sides of the first jet and are called the second jet. Then, the third pulsation results in the generation of the third jets between the first and second jets. The first jet is similar to a highly focused supersonic microjet in other studies (Tagawa et al. 2012;Peters et al. 2013). Tagawa et al. (2012) argued that the supersonic microjet is induced by 'kinematic focusing' of the liquid. This is because the impulsive acceleration or pressure converges toward the centre of curvature of the droplet surface. A shock wave is the dominant factor in the kinematic focusing of the liquid when it is reflected at the surface in their study. However, in the present study, when the generated shock wave reaches the droplet surface, the pressure behind it is of the order of dozens of MPa, even several MPa according to the study of Wang et al. (2018b). Furthermore, the droplet surface can achieve great acceleration owing to the initial expansion of cavitation bubble, as shown in figure 10(b). Consequently, the bubble expansion may be one of the most important factors in generating the first jet in the present study. Figure 16(a) shows the velocity field of a droplet with a diameter of 12 mm during the first pulsation of the cavitation bubble obtained with numerical simulation. As shown in figure 16(a-2), we observe the focusing of the flow at the droplet as marked by the red  Figure 16. Evolution of the ring-shaped jet for a droplet with a diameter of 12 mm obtained numerically: (a) velocity field; (b) velocity at the surface centre. In (a), the cross section is at z = 0.4 mm, and the focusing of the flow at the droplet surface is induced by the propagation of the shock wave, as indicated by the red ellipse. In (b), the velocity at the centre of the droplet surface is that of the ring-shaped jet. ellipse, which can be induced by the shock wave according to its propagation velocity. This indicates that the initial perturbation is stimulated by the shock wave. As the bubble expands, the droplet surface grows gradually, forming the first jet. Given the velocity curve of the jet in figure 16(b), the droplet surface is accelerated during two periods. From the perspective of the time scale, it is found that Period 1 is mainly caused by the propagation of shock wave, while Period 2 is dominated by the expansion of the bubble. Subsequently, Period 3 shows that the droplet surface decelerates owing to the surface tension until a small droplet is generated at the tip of the first jet at approximately 600 μs. As a result, the generation and acceleration of the jet is caused by the combination of the shock wave and expansion of the cavitation bubble.