1. Introduction
When a liquid interface undergoes periodic oscillations with sufficient amplitude, standing-wave patterns are generated. This phenomenon is known as Faraday instability, named after its discoverer Faraday (Reference Faraday1831), who first observed its occurrence on the free surface of a vertically vibrating liquid bath. Faraday also found that the resulting wave patterns oscillate at half the frequency of the driving force. This type of response is denoted as half-harmonic, and is a hallmark of this kind of instability. Faraday instability is also referred to as parametric instability because it involves the variation of a key parameter in the dynamical system – in this case, the effective gravitational acceleration. Faraday’s findings were later validated by Rayleigh (Reference Rayleigh1883a , Reference Rayleighb ). The first linear stability analysis of Faraday waves for an inviscid liquid with a flat surface was performed by Benjamin & Ursell (Reference Benjamin and Ursell1954), who demonstrated that the wave evolution follows a Mathieu equation. This stability analysis shows that wave frequencies can be half-harmonic, harmonic, ultra-harmonic (integer multiples of half the driving frequency) or super-harmonic (integer multiples of the driving frequency). However, the extension of this analysis to viscous fluids by Kumar & Tuckerman (Reference Kumar and Tuckerman1994) revealed that the onset of the half-harmonic response requires the lowest acceleration, making it the only response that manifests. Weak nonlinear effects can lead to interactions between different surface-wave modes, influencing the selection of the dominant wave pattern. Fauve (Reference Fauve1998) demonstrated that in an axisymmetric vessel, the most unstable linear eigenmode is an axisymmetric standing wave. However, this axisymmetric pattern is nonlinearly unstable, and nonlinear effects ultimately favour the emergence of non-axisymmetric patterns. These patterns can include lattices of stripes, squares or hexagons (Kudrolli & Gollub Reference Kudrolli and Gollub1996), as well as more complex structures such as quasi-patterns (Christiansen, Alstrom & Levinsen Reference Christiansen, Alstrom and Levinsen1992; Edwards & Fauve Reference Edwards and Fauve1994), triangles (Müller Reference Müller1993), superlattice patterns (Kudrolli & Gollub Reference Kudrolli and Gollub1996; Wagner, Müller & Knorr Reference Wagner, Müller and Knorr1999), oscillons (Arbell & Fineberg Reference Arbell and Fineberg2000a ) and rhomboidal patterns (Arbell & Fineberg Reference Arbell and Fineberg2000b ). The influence of surfactants on the interface of liquids with arbitrary depth was investigated by Kumar & Matar (Reference Kumar and Matar2002b , Reference Kumar and Matar2004) and extended to the specific case of thin liquid layers by Kumar & Matar (Reference Kumar and Matar2002a ) and Matar, Kumar & Craster (Reference Matar, Kumar and Craster2004).
Faraday instability on spherical fluid interfaces was first examined by Kelvin (Reference Kelvin1863) and later by Rayleigh (Reference Rayleigh1879). A key feature of the Faraday instability is the direct relationship between a specific pattern wavenumber and its frequency, as described by a dispersion relation. In the case of a spherical interface, this relationship was derived by Rayleigh (Reference Rayleigh1879) and Lamb (Reference Lamb1932), by considering the restoring force due to surface tension in a system of two incompressible, inviscid fluids, where one fluid entirely envelops the other, and reads

where
$\omega _{0,l}$
is the angular frequency of the pattern,
$l$
is its angular wavenumber,
$\rho _{ { i}}$
and
$\rho _{ { e}}$
are the densities of the internal and external fluids, respectively,
$\sigma$
represents surface tension and
$R_0$
is the radius of the spherical interface between the two fluids. Due to the spherical symmetry, the wavenumber must take on integer values. Consequently, when the driving frequency
$\omega _{ { d}} = 2\omega _{0,l}$
corresponds to a non-integer wavenumber according to the dispersion relation (1.1), the actual wavenumber adjusts to the nearest integer, leading to a quantised response spectrum. However, in such cases, the onset acceleration required to generate a pattern will exceed the value necessary when the driving frequency directly resonates with an integer wavenumber.
The interface deformation,
$\eta (\theta , \phi , t)$
, is classically expressed using real spherical harmonics,
$Y_{l}^{m}(\theta , \phi )$
, as follows:


where
$\theta$
and
$\phi$
are the polar and azimuthal angles in spherical coordinates, respectively;
$P_{l}^{m}(\cos (\theta ))$
denotes the associated Legendre function of degree
$l$
and order
$m$
;
$N_{l}^{m} = (\max (P_{l}^{m} (\cos (\theta )) \cos (m \phi )))^{-1}$
is a normalisation factor that ensures the value of the spherical harmonic is constrained to unity; and
$a_l^{m}(t)$
represents the time-dependent amplitude of each spherical harmonic. By normalising the harmonics in this way, the shape mode magnitudes,
$a_l^m(t) Y_l^m$
, can be directly compared through the coefficients
$a_l^m(t)$
alone. Spherical harmonics are also known as shape modes and are ordered by their wavenumber pairs [
$l, m$
]. For a given
$l \neq 0$
, these modes are classified as follows: zonal modes when
$m = 0$
, sectoral modes when
$m = l$
and tesseral modes when
$m \neq 0$
and
$m \neq l$
. Figure 1, which illustrates modes up to the sixth degree, shows that zonal modes are axisymmetric with
$l$
azimuthal nodal lines, sectoral modes exhibit star symmetry with
$m$
polar nodal lines and tesseral modes have (
$l - m$
) azimuthal nodal lines along with
$m$
polar nodal lines. The
$l=0$
mode represents a uniform spherical deformation, often referred to as breathing mode. It characterises radial oscillations of bubbles driven by a sinusoidal acoustic field. The
$l=1$
mode manifests as an alternating translational motion and must correspond to a gravity wave rather than a capillary wave, as it does not deform the interface and thus cannot be restored by surface tension. It is crucial to note that the dispersion relation (1.1), which describes the linear stability of the interface, is independent of the spherical harmonic order
$m$
. Consequently, a departure from the spherical symmetry state leads to a set of
$l+1$
linearly independent and equally probable solutions. In this context, the spectrum is termed as degenerate. The final deformation pattern, determined by a specific combination of shape modes with the same
$l$
as represented in (1.2), will be selected by nonlinear effects (Chossat, Lauterbach & Melbourne Reference Chossat, Lauterbach and Melbourne1991).

Figure 1. Shape modes of a spherical interface ordered by degree
$l$
and order
$m$
, up to the sixth degree. Zonal modes occupy the column marked in red (
$m=0$
), sectoral modes are located along the highlighted diagonal in blue (
$l=m$
) and tesseral modes are distributed across the remaining yellow area (
$l\neq m$
). The
$l=0$
mode corresponds to a uniform spherical deformation and is commonly referred to as breathing mode.
Drops and bubbles are two of the most common forms of fluids encapsulated within one another. The predictions of the Rayleigh–Lamb spectrum for small shape oscillations of liquid drops have been experimentally confirmed by Trinh, Zwern & Wang (Reference Trinh, Zwern and Wang1982) by suspending the drops in a host fluid through acoustic forces. To eliminate the influence of acoustic trapping and enable fully non-invasive investigations, researchers have also conducted experiments in space using free-floating drops (Wang, Anilkumar & Lee Reference Wang, Anilkumar and Lee1996) or during parabolic flights employing a liquid layer enclosed within a spherical container (Falcón et al. Reference Falcón, Falcon, Bortolozzo and Fauve2009). Nevertheless, the applicability of the Rayleigh–Lamb spectrum is limited to inviscid fluids and infinitesimal interfacial deformations. The influence of viscosity has been studied by Miller & Scriven (Reference Miller and Scriven1968), Prosperetti (Reference Prosperetti1980) and Ebo Adou & Tuckerman (Reference Ebo Adou and Tuckerman2016), who found that increased viscosity raises the acceleration threshold for the onset of shape modes, particularly as the wavenumber increases. Concerning finite deformations, Tsamopoulos & Brown (Reference Tsamopoulos and Brown1983) were the first to analytically investigate weakly nonlinear effects on inviscid axisymmetric drops. Later, Lundgrent & Mansour (Reference Lundgrent and Mansour1988) employed a boundary integral method to compute large axisymmetric motions with weak viscous effects. More recently, Ebo-Adou et al. (Reference Ebo-Adou, Tuckerman, Shin, Chergui and Juric2019) conducted fully three-dimensional viscous simulations, which revealed the final deformation patterns of the interface as determined by nonlinear effects. The
$l = 2$
mode alternates between oblate and prolate spheroidal shapes, and modes with
$l \geqslant 3$
produce patterns analogous to Platonic solids. Building on this, Panda et al. (Reference Panda, Kahouadji, Abdal, Tuckerman, Shin, Chergui, Juric and Matar2024) extended the study to higher spherical harmonic degrees and amplitudes.
The Rayleigh–Lamb spectrum, originally developed for free spherical interfaces, has been extended to also account for constrained spherical inviscid drops. Strani & Sabetta (Reference Strani and Sabetta1984) first introduced this extension by considering a constraint in the form of a spherical bowl. Later, Bostwick & Steen (Reference Bostwick and Steen2014) explored the dynamics of an inviscid sessile drop resting on a flat, planar support, analysing its natural oscillations. They showed that these oscillations are influenced by the static contact angle and the mobility of the contact line. Notably, they also found that a flat support breaks the spectral degeneracy, leading to shape modes with the same degree
$l$
but different order
$m$
exhibiting distinct resonance frequencies. These theoretical predictions align reasonably well with experimental observations by Chang et al. (Reference Chang, Bostwick, Steen and Daniel2013, Reference Chang, Bostwick, Daniel and Steen2015) for pinned drops with minimal gravitational effects, although the agreement weakens for flatter drops and higher-order modes. The researchers also discovered the phenomenon of mode mixing, where multiple mode shapes are excited by a single driving frequency. An alternative analytical approach for predicting the shape-mode spectrum of sessile drops, proposed by Sharma & Wilson (Reference Sharma and Wilson2021), provides more accurate results for flatter drops. In their experimental study, Vukasinovic, Smith & Glezer (Reference Vukasinovic, Smith and Glezer2007) found that harmonic axisymmetric standing waves emerge even under minimal forcing. When the forcing amplitude exceeds a critical threshold, the Faraday instability generates half-harmonic azimuthal waves along the contact line. At even higher forcing amplitudes, these waves merge with the harmonic waves, creating a lattice-like pattern. These distinctive patterns were also reported in three-dimensional viscous simulations by Panda et al. (Reference Panda, Kahouadji, Tuckerman, Shin, Chergui, Juric and Matar2023).
Compared with drops, bubbles exhibit richer dynamics owing to the additional degree of freedom introduced by the compressibility of the inner fluid. When subjected to an acoustic driving with a wavelength significantly larger than the bubble size, free bubbles undergo spherical oscillations, which correspond to the
$l=0$
mode in spherical harmonics terminology. The natural angular frequency
$\omega _{0,0}$
of this mode, for small deformations, can be determined by linearising the Rayleigh–Plesset equation in combination with the polytropic gas approximation, which describes the dynamics of a spherical bubble (Plesset & Prosperetti Reference Plesset and Prosperetti1977), as follows:

where
$\rho _l$
is the density of the liquid medium,
$n$
is the gas polytropic index and
$p_{\infty }$
is the ambient pressure. In addition to shape oscillations induced by direct forcing, as seen with drops, bubbles can experience shape oscillations induced by the harmonic radial acceleration stemming from their spherical oscillations. To predict the amplitude of shape oscillations in the context of inviscid fluids, a classical approach involves employing a second-order ordinary differential equation linearised with respect to the mode amplitude (Plesset Reference Plesset1954; Benjamin & Strasberg Reference Benjamin and Strasberg1958; Hsieh & Plesset Reference Hsieh and Plesset1961; Eller & Crum Reference Eller and Crum1970). The corresponding problem for viscous fluids has been treated by Prosperetti (Reference Prosperetti1977), Ceschia & Nabergoj (Reference Ceschia and Nabergoj1978), Brenner, Lohse & Dupont (Reference Brenner, Lohse and Dupont1995) and Hilgenfeldt, Lohse & Brenner (Reference Hilgenfeldt, Lohse and Brenner1996). By performing a first-order analysis, Francescutto & Nabergoj (Reference Francescutto and Nabergoj1978) derived a threshold condition for the onset of shape modes on an oscillating bubble in a sound field. Versluis et al. (Reference Versluis, Goertz, Palanchon, Heitman, Van Der Meer, Dollet, De Jong and Lohse2010) leveraged these analytical findings to compare the theoretical acoustic pressure thresholds for the onset of shape modes against experimental results obtained from testing free micrometric bubbles driven by ultrasound, finding a good agreement.
However, these studies do not account for the nonlinear impact of a specific shape mode on the breathing mode, bubble translation and other shape modes. Ignoring this influence in the aforementioned models results in predictions of unbounded growth for any shape deformation. As a result, the applicability of these models is limited to predicting the amplitude threshold for the onset of shape oscillations. Initially, the coupling between the breathing mode and a single shape mode was studied analytically by Longuet-Higgins (Reference Longuet-Higgins1989a , Reference Longuet-Higginsb , Reference Longuet-Higgins1991), Mei & Zhou (Reference Mei and Zhou1991) and Feng & Leal (Reference Feng and Leal1992, Reference Feng and Leal1994, Reference Feng and Leal1997) and shown to be both periodic and energy conserving. Building on this, subsequent research has progressively advanced the understanding of the interaction between the breathing mode, translational motion and multiple shape modes of a bubble (Feng & Leal Reference Feng and Leal1995; Doinikov Reference Doinikov2004; Shaw Reference Shaw2006, Reference Shaw2009, Reference Shaw2017). Experimental studies with acoustically trapped micrometric bubbles driven by ultrasound carried out by Guédra et al. (Reference Guédra, Inserra, Mauger and Gilles2016, Reference Guédra, Cleve, Mauger, Blanc-Benon and Inserra2017) have confirmed these theories, reporting strong nonlinear mode coupling, including the saturation of instability and the triggering of non-parametric shape modes. The development of nonlinear coupled models has been pivotal in elucidating the mechanisms behind the sound emission stemming from the activation of the breathing mode (Minnaert Reference Minnaert1933) and the intriguing ‘dancing motion’ of a bubble in an acoustic field (Kornfeld & Suvorov Reference Kornfeld and Suvorov1944). Cattaneo et al. (Reference Cattaneo, Guerriero, Shakya, Krattiger, G. Paganella, Narciso and Supponen2025) experimentally characterised the pattern for the first six shape modes of ultrasound-driven phospholipid-coated microbubbles in contact with a supersoft, superhydrophilic substrate. This configuration minimises the effects of wall confinement, allowing the bubbles to behave as if they are in an unbounded fluid. The shape-mode patterns found correspond to the numerical predictions from Ebo-Adou et al. (Reference Ebo-Adou, Tuckerman, Shin, Chergui and Juric2019) for droplets, suggesting a common interfacial response to both drops and bubbles, at least when the two cases involve similar density and viscosity ratios, as well as comparable surface tension.
While we have accumulated substantial knowledge regarding the oscillations of bubbles in unbounded fluids, our understanding of how confining surfaces influence bubble dynamics remains relatively limited. Strasberg (Reference Strasberg1953) was the first to quantify the influence of an infinitely rigid wall on the breathing mode resonance of a single bubble using the method of images. Shklyaev & Straube (Reference Shklyaev and Straube2008) and Fayzrakhmanova, Straube & Shklyaev (Reference Fayzrakhmanova, Straube and Shklyaev2011) explored the theoretical aspects of linear natural and forced shape oscillations of a hemispherical bubble on a solid substrate, examining the effects of bubble compressibility and contact angle hysteresis. Their findings revealed that the contact line dynamics induces an interaction between the breathing mode and the shape modes, even within a linear framework. Prosperetti (Reference Prosperetti2012) and Vejrazka, Vobecka & Tihon (Reference Vejrazka, Vobecka and Tihon2013) investigated the shape modes of a spherical bubble pinned to a ring, while Ding & Bostwick (Reference Ding and Bostwick2022) analysed the hydrodynamic stability of a bubble on a solid substrate. That study uncovered a complex relationship between the frequency spectrum and wetting properties, which is influenced by factors such as the static contact angle, the equilibrium bubble pressure and the contact-line dynamics, and showed that the constraint removes the mode degeneracy seen in free bubbles. Finally, Fauconnier, Béra & Inserra (Reference Fauconnier, Béra and Inserra2020) investigated experimentally the modal behaviour of a wall-attached bubble under acoustic excitation. Through the analysis of top-view recordings using spherical harmonics decomposition, the study suggested the absence of mode degeneracy, as well as the presence of modal mixing and competition.

Figure 2. (a) Sessile water droplet with a 5 mm radius subjected to vertical vibrations at a frequency of 1080 Hz. Reproduced with permission from the work by Vukasinovic, Glezer & Smith (Reference Vukasinovic, Glezer and Smith2000). (b) Wall-attached air bubble with a 68
$\;\unicode{x03BC}$
m radius subjected to ultrasound driving at a frequency of 100 kHz. Reproduced with permission from the work by Cattaneo et al. (Reference Cattaneo, Presse, Shakya and Supponen2023).
When the interface acceleration surpasses a certain threshold, the amplitude of Faraday waves rapidly increases, causing the wave profile to become unstable. This instability ultimately leads to the formation of a jet that moves from the higher-density fluid to the lower-density fluid. Longuet-Higgins (Reference Longuet-Higgins1983) was the first to observe this phenomenon on a vertically vibrating flat water–air surface. In some remarkable cases, he observed jets produced by a relatively mild excitation amplitude of 0.5 mm that rose up to a height of 1.7 m. The understanding of Faraday wave-induced jets on flat interfaces was significantly advanced by the work of Goodridge et al. (Reference Goodridge, Shi and Lathrop1996, Reference Goodridge, Tao, Hentschel, H.G. and Lathrop1997), Tao Shi, Goodridge & Lathrop (Reference Tao Shi, Goodridge and Lathrop1997), Hogrefe et al. (Reference Hogrefe, Peffley, Goodridge, Shi, Hentschel and Lathrop1998), Goodridge, Hentschel & Lathrop (Reference Goodridge, Hentschel and Lathrop1999) and Zeff et al. (Reference Zeff, Kleber, Fineberg and Lathrop2000). These studies identified the critical wave height necessary for jet formation and demonstrated that the surface undergoes a self-similar collapse prior to jetting. Vukasinovic et al. (Reference Vukasinovic, Glezer and Smith2000, Reference Vukasinovic, Glezer and Smith2001) and James et al. (Reference James, Vukasinovic, Smith and Glezer2003) were the first to observe Faraday wave-induced jetting in a spherical sessile water drop placed on a vibrating plate. A sequence of images from their studies on drop jetting is reproduced in figure 2(a). The jets, which extended outward from the droplet, led to atomisation and drop bursting through Rayleigh–Plateau instability. Recently, Prasanna et al. (Reference Prasanna, Biasiori-Poulanges, Yu, El-Rabii, Lukić and Supponen2024) studied viscoelastic drops and demonstrated that the high resistive stresses in the jets significantly enhance their stability against the Rayleigh–Plateau instability. Cattaneo et al. (Reference Cattaneo, Presse, Shakya and Supponen2023) provided visualisations that clearly depict the formation of cyclic jets fostered by the Faraday instability on wall-attached bubbles oscillating under ultrasound driving. Figure 2(b) presents a series of images capturing this process. We believe that cyclic jets driven by Faraday instabilities, though previously unrecognised, have likely been observed in multiple earlier studies employing driving frequencies spanning from hertz (Crum Reference Crum1979), through kilohertz (Prabowo & Ohl Reference Prabowo and Ohl2011), and up to megahertz (Vos et al. Reference Vos, Dollet, Versluis and De Jong2011). Shape-mode-induced jets were also observed for bubbles in contact with a vibrating surface, as reported by Biasiori-Poulanges et al. (Reference Biasiori-Poulanges, Bourquard, Lukić, Broche and Supponen2023). Recently, Dhote et al. (Reference Dhote, Kumar, Kayal, Goswami and Dasgupta2024) successfully used Gerris, an open-source computational fluid dynamics software, to numerically simulate and demonstrate the formation of a shape-mode-induced jet in an inviscid, incompressible, pre-deformed bubble. Finally, in our previous work (Cattaneo et al. Reference Cattaneo, Guerriero, Shakya, Krattiger, G. Paganella, Narciso and Supponen2025), we studied the jet formation driven by Faraday instability on practically unconstrained, phospholipid-coated, microbubble ultrasound contrast agents subjected to high-frequency ultrasound, revealing a multi-jetting behaviour, where the number of jets scales with the number of lobes in the shape mode. We also discovered that these jets are capable of piercing cell membranes, enabling precise and targeted drug delivery.
While significant progress has been made in understanding oscillations and jetting in acoustically driven free bubbles, the dynamics of wall-attached bubbles remains poorly understood. In particular, the selection and evolution of shape-mode patterns in these bubbles have not been experimentally resolved in three dimensions, and the resulting shape-mode-induced jetting has yet to be compared with that observed in free bubbles. The present study addresses this gap by examining bubbles, with radii of around 100
$\;\unicode{x03BC}$
m, in contact with a rigid substrate, subjected to low-frequency ultrasound. Building on the work of Fauconnier et al. (Reference Fauconnier, Béra and Inserra2020), which was limited to top-view imaging, we introduce a dual-imaging technique that combines bright-field microscopy with phase-contrast X-ray synchrotron imaging. This approach captures bubble dynamics from two orthogonal perspectives, providing visual access to the full three-dimensional bubble shape. The X-ray modality minimises refraction artefacts and reveals all bubble folds, including those on the distal side that are hidden in traditional back illumination. This dual-view methodology enables, for the first time, direct observation of the formation, spectrum and temporal evolution of Faraday instability leading up to jet formation in wall-attached bubbles. The results allow quantitative comparison with theoretical predictions and three-dimensional numerical simulations, offering new insights into the mechanisms driving these complex phenomena. This paper is structured as follows. Section 2 outlines our experimental set-up that combines X-ray synchrotron and conventional light sources to enable dual-view imaging. Section 3 outlines the theoretical framework for describing bubble deformations in contact with a wall and explores the permissible shape modes as a function of the wetting conditions. Section 4 provides an overview of the bubble response, breaking it down into four sequential regimes. Section 5 presents a quantitative experimental analysis of the four regimes, supplemented by comparisons with theoretical predictions and numerical simulations, and offers physical interpretations of the findings. Section 6 provides a comparison between the experimental visualisations and three-dimensional simulations of bubble dynamics generated using the boundary element method (BEM).
2. Experimental set-up
A schematic of the experimental set-up is shown in figure 3(a). In a water bath
$(T_l\approx {22}\,^\circ {\textrm{C}} )$
filled with deionised water, a rising stream of monodisperse bubbles with an equilibrium radius within the range
$R_0={60}{-}{140}\;\unicode{x03BC}{\textrm{m}}$
is generated using a polydimethylsiloxane microfluidic chip fabricated using soft lithography. The chip layout consists of a T-junction where the channel conveying the continuous phase (water) intersects with the channel transporting the dispersed phase (air), resulting in the formation of bubbles. Additionally, two channels further downstream create a sheath flow which assist the bubble separation by increasing the flow rate. The geometry of the microfluidic chip is depicted in figure 3(b). The liquid phase is supplied by a syringe pump (NE-300, Darwin Microfluidics) with a flow rate of
${50} \;\unicode{x03BC} {\textrm {l min}}^{-1}$
. The gas is fed by an air compressor (Fatmax DST 101/8/6, Stanley) and its pressure, which determines the bubble size, is controlled with a pressure-reducing valve (RP1000-8G-02, CKD). Finally, the sheath flow is provided by a syringe pump (NE-300, Darwin Microfluidics) with a combined flow rate of 200
$\;\unicode{x03BC} {\textrm {l min}}^{-1}$
.

Figure 3. (a) Schematic of the experimental set-up. (C1, C2) Cameras, (GC) glass capillary, (LI) LED illuminator, (MC) microfluidic chip, (OL1, OL2) objective lenses, (S) scintillator, (SH) sample holder, (SR) sound reflector, (TL1, TL2) tube lenses, (TW1, TW2) telescopic windows, (US) ultrasound transducer. (b) Geometry of the microfluidic bubble-generator chip. A single bubble is diverted from the upward bubble stream and propelled towards the bottom of the glass capillary using a manually operated syringe. (c) Testing conditions with a single microbubble positioned underneath a glass capillary. (d) Acoustic driving pressure produced by the ultrasound transducer and measured by a needle hydrophone, normalised to the steady-state amplitude value.
One of these micrometric air bubbles is diverted from the stream using a manually operated syringe and placed to rest on the bottom of a square borosilicate glass capillary immersed in water. A square capillary, characterised by a side length
$l$
of 1.4 mm and a wall thickness
$\delta$
measuring 0.2 mm (8100-50, CM Scientific), is chosen over a solid substrate because it mitigates sound reflections due to its hollow structure filled with water and capillary walls significantly thinner than the wavelength
$\lambda$
of the acoustic driving (
$ \lambda / \delta = 250$
). Nevertheless, its side length is sufficiently large to avoid boundary effects (
$l/ 2R_0 \gt 5$
). Moreover, the capillary cross-section provides good bending rigidity to the substrate, preventing unintended oscillations when subjected to acoustic stimulation. Figure 3(c) presents a close-up view of the bubble in contact with the capillary. The bubble is acoustically driven at a frequency
$f_{ { d}} = {30}\,{\textrm{kHz}}$
with pressure amplitudes in the range
$p_{{ a}}={1}{-}{15} \,{\textrm{kPa}}$
using an ultrasound transducer (GS30-D25, The Ultran Group) positioned in the water bath perpendicular to the horizontal plane. The driving pulse is generated using a function generator (LW420B, Teledyne LeCroy) and subsequently amplified by a radiofrequency power amplifier (1020L, E&I). The resulting ultrasound pulse shape, as measured by a needle hydrophone (
${0.2}\, {\textrm{mm}}$
, NH0200, Precision Acoustics), is depicted in figure 3(d).
Side-view visualisations of the bubble response are acquired using high-speed synchrotron X-ray phase-contrast radiography at the ID19 beamline of the European Synchrotron Radiation Facility (ESRF), during the
$7/8 +1$
filling mode with an integral storage ring current of
${200}\;\textrm {mA}$
. We employ X-ray synchrotron radiation because it is far less susceptible to refraction on curved surfaces compared with visible light. This choice helps eliminate the black bands typically seen on spherical surfaces illuminated with visible light (as observed in figure 2
b). These bands can obscure the surface evolution of the bubble during lobe folding and hinder the early detection of jet formation. Moreover, the quasi-parallel X-ray beam at the ID19 beamline, enabled by its exceptionally long source-to-sample distance, effectively minimises source size effects (penumbral blur), enabling us to clearly resolve all interface folds, even on the distal side of the bubble. This enhanced clarity makes it easier to identify the shape modes that develop, particularly when the bubble loses its axisymmetry. X-ray light requires the use of indirect detectors that first convert X-rays into visible light and then into electrons. However, these detectors lose efficiency as they approach the optical resolution limit, necessitating ultrathin scintillators for effective performance. As a result, high-speed X-ray imaging often comprises spatial resolution to enhance the signal-to-noise ratio, leading to less sharply defined interfaces compared with those commonly seen in visible-light imaging.
For this study, the bubble in contact with the glass capillary is illuminated by a pink X-ray beam generated by the U17.6 undulator of the beamline, which is set with a 12 mm gap (undulator period
$\lambda _{ { x}} = {17.6}\,\textrm {mm}$
, number of periods of magnet pairs
$N_{ { x}} = 92$
). The insertion device produces an on-axis central beam energy (first harmonic) at 17.9 keV, with a
$2\,\%$
bandwidth (full width at half maximum), while the second and third harmonics are suppressed by two orders of magnitude in photon flux. The beam is filtered with a set of mandatory optical elements along the vacuum flight path: a 1.8 mm thick diamond window, a 0.7 mm thick aluminium filter and successive thin carbon and beryllium windows. The beam is collimated to the field of view by using eight compound refractive lenses from the beamline transfocator, located approximately 35 m from the source, acting as beam condensers. Additionally, two sets of in-vacuum slits are employed to crop the beam to match the field of view. The indirect detector consists of a 250
$\;\unicode{x03BC}$
m thick LuAG:Ce scintillator (Crytur), coupled via a mirror to a microscope objective lens of focal length
${\mathcal{F}} = {10} \,\textrm {mm}$
(MY20X-804, Mitutoyo) with lead glass protection and a
${\mathcal{F}} = {200}\, {\textrm{mm}}$
tube lens for a total magnification of
${\times}20$
. The luminescent image is then relayed to a high-speed camera (HPV-X2, Shimadzu). The distance between the sample and the scintillator is around 3 m.
The water tank is constructed from polymethylmethacrylate (PMMA) and features two water-tight telescopic windows on two opposite-facing walls comprising aluminium tubes (
$\varnothing 0.5^{\prime\prime}$
stackable tube lens, Thorlabs) and circular 0.5 mm thick PMMA windows secured at their ends by retaining rings and sealed with high-vacuum grease. These telescopic windows are designed to minimise the path traversed by X-ray radiation in water, aiming to decrease absorption and consequently enhance the signal. Therefore, the distance between the windows is minimised to the greatest extent, compatible with spatial limitations dictated by other components, at approximately
$d= {10}\,\textrm {mm}$
. The transmission coefficient of the X-ray signal, pertaining to the thickness of the water layer utilised, can be estimated to be approximately equal to
$T=0.4$
using the relation
$T = {\textrm e}^{- \epsilon _{ { l}} \rho _{ { l}} d}$
(Henke, Gullikson & Davis Reference Henke, Gullikson and Davis1993), where
$\epsilon _{ { l}}$
is the mass absorption coefficient, which for water at
${17.9}\,\textrm {keV}$
is measured to be
${9.263}\times {10^{-2}}\,{\textrm {m}^2}\,{\textrm {kg}^{-1}}$
(Hubbell & Seltzer Reference Hubbell and Seltzer1995), and
$\rho _{ { l}}$
is the water density.
Simultaneous top-view visualisations of the bubble response are captured using high-speed visible-light microscopy. A custom-built upright microscope is realised for the purpose using modular optomechanics components (Thorlabs, cage system) and equipped with a water-dipping objective lens of focal length
${\mathcal{F}} = {20} \,\textrm {mm}$
(N10XW-PF, Nikon) and a
${\mathcal{F}} = {400} \,\textrm {mm}$
tube lens (TL400-A, Thorlabs) for a total magnification of
${\times}20$
. Backlight illumination is provided by a custom-built high-power LED illuminator used in continuous mode for live imaging and in burst mode for video recording. The light emission is guided with a liquid optical fibre and sent through the sample and then into the objective by reflecting it against the upper surface of the ultrasound transducer. Sound reflections against the objective lens are mitigated by attaching a transparent PMMA prism to its base (see figure 3
c). Video recordings are captured using a high-speed camera (HPV-X2, Shimadzu).
The top-view and side-view cameras are synchronised (
${\lt}{10}\,\textrm {ns}$
of delay between the two) and allow for recordings at a frame rate of 150 kHz over 256 frames of continuous visualisation of a
${640\;\unicode{x03BC}\textrm {m} \times 400}\;\unicode{x03BC}{\textrm{m}}$
field of view with a 1.6
$\;\unicode{x03BC}$
m pixel resolution. The activation of the ultrasound pulse, camera recording, light flash and X-ray beam switching are synchronised through a delay generator (DG645, Stanford Research Systems).
Due to X-rays having a shorter wavelength than visible light, radiographs exhibit a remarkably high depth of field. While this ensures sharp focus across all image planes, it impedes gaugeing the distance of an object from the optical window through distance blur. Consequently, positioning the bubble under study within the acoustic focal volume requires the combined use of side and top views. The one-time alignment protocol implemented to ensure the coincidence of the centres of the fields of view and the acoustic focal point is as follows. First, we align the X-ray optical path within the telescopic windows by repositioning the optical table supporting the set-up. Next, we position the needle hydrophone parallel to the horizontal plane passing through the X-ray optical path, halfway between the two telescopic windows, and align its tip with the centre of the X-ray field of view through a motorised three-axis microtranslation stage (PT3/M-Z8, Thorlabs). We manoeuvre the ultrasound transducer using a manual three-axis microtranslation stage (three DTS50/M, Thorlabs) to obtain the maximum signal from the hydrophone. We position the entire top-view videomicroscopy system by means of a manual three-axis microtranslation stage (12694, 66-511, Edmund Optics) so that the tip of the hydrophone is in the centre of the microscope field of view and in focus. Finally, the hydrophone is replaced by the square capillary.

Figure 4. Definition sketch of a bubble resting against a rigid flat substrate with a static contact angle
$\alpha _0$
. The cross-sections along the planes defined by projection lines (a)
$\text{A}{-}\text{A}$
and (b)
$\text{B}{-}\text{B}$
. The equilibrium surface, represented by a dashed line, is defined as
${\boldsymbol{R}}_0(\theta ,\phi )$
, where polar angle
$\theta \in [-\alpha _0,\alpha _0]$
and the azimuthal angle
$\phi \in [0, 2\pi ]$
serve as surface coordinates. The time-dependent surface deformation is denoted as
$\delta {\boldsymbol{R}}(\theta ,\phi ,t)$
. The deformed surface
${\boldsymbol{R}}_0(\theta ,\phi ) + \delta {\boldsymbol{R}}(\theta ,\phi ,t)$
is depicted by a solid line. Vectors
$\boldsymbol{n}$
,
$\boldsymbol{t}$
and
$\boldsymbol{b}$
are the normal, tangential and binormal unit vectors to the equilibrium surface, respectively. The contact line
$\boldsymbol{\Gamma }$
is shown as a dotted line. The surface deformation normal to the equilibrium surface
$\delta {\boldsymbol{R}}(\theta ,\phi ,t) \,{\boldsymbol\cdot}\, \boldsymbol{n} = \eta$
in the (c)
$\text{A}{-}\text{A}$
and (d)
$\text{B}{-}\text{B}$
planes.
3. Theoretical framework
3.1. Kinematic model
Let us consider a bubble resting against a rigid flat substrate, as depicted in the side and top views in figure 4(a,b). Given the small Bond number,
$Bo = {(\rho _{ { l}} - \rho _{ { g}}) g R_0^2}/{\sigma } \sim O(10^{-3})$
, where
$(\rho _{ { l}} - \rho _{ { g}})$
is the density difference between water and air and
$g$
is the gravitational acceleration, the influence of gravitational effects is negligible for the problem at hand. Consequently, under equilibrium conditions, the bubble can be approximated as a spherical cap with a static contact angle
$\alpha _0$
. We define the origin of the reference frame at the geometric centre of the uncut sphere. The surface of the bubble at equilibrium can be parametrically described as
${\boldsymbol{R}}_0(\theta ,\phi )$
, where
$\theta \in [-\alpha _0,\alpha _0]$
is the polar angle and
$\phi\,{\in}\,[0, 2\pi ]$
is the azimuthal angle. By normalising the bubble equilibrium surface (
$\|{\boldsymbol{R}}_0\| = 1$
), the arc angle and arc length values are equal. The surface deformation is defined as
$\delta {\boldsymbol{R}}(\theta ,\phi ,t)$
. This deformation may displace the contact line, denoted as
$\boldsymbol{\Gamma }(\theta )$
, along the substrate. The unit vector normal to
${\boldsymbol{R}}_0$
is denoted by
$\boldsymbol{n}$
. The unit vector tangential to
${\boldsymbol{R}}_0$
and normal to
$\boldsymbol{\Gamma }$
is denoted by
$\boldsymbol{t}$
. The unit vector tangential to
${\boldsymbol{R}}_0$
and to
$\boldsymbol{\Gamma }$
is denoted by
${\boldsymbol{b}} = {\boldsymbol{n}} \times {\boldsymbol{t}}$
. The unit vector normal to the substrate is denoted as
${\boldsymbol{n}}_{\textrm {s}}$
. The bubble deformation can be decomposed into the component normal to the bubble equilibrium surface, defined as

and the component tangential to the bubble equilibrium surface, represented by

as
$\delta {\boldsymbol{R}} \perp {\boldsymbol{b}}$
. An example of normal deformation
$\eta$
in the plane perpendicular to the substrate defined by the projection line
$\text{A}{-}\text{A}$
, and thus identified as
$\eta (\theta ,0,t)$
, is shown in figure 4(c). In contrast, an example of normal deformation
$\eta$
in the plane parallel to the substrate defined by the projection line
$\text{B}{-}\text{B}$
, and thus identified as
$\eta (\pi /2,\phi ,t)$
, is shown in figure 4(d).
The equilibrium contact angle
$\alpha _0$
is related to the surface tensions of the solid–liquid
$\sigma _{{ sl}}$
, solid–gas
$\sigma _{\textit{sg}}$
and liquid–gas
$\sigma _{{ lg}}$
through the Young–Dupré equation:

A perturbation in the bubble surface
$\delta {\boldsymbol{R}}$
results in a change in the wetting condition (Tyuptsov Reference Tyuptsov1966; Myshkis et al. Reference Myshkis, Babskii, Kopachevskii, Slobozhanin and Tyuptsov1987):

A full derivation of this relation is provided in Appendix A. We use this kinematic boundary condition to determine the spectrum of shape modes for the deformation
$\delta {\boldsymbol{R}}$
under various wetting conditions.
3.2. Shape modes for wall-attached bubbles
The tangential deformation
$\tau (\theta ,\phi ,t)$
does not affect the curvature of the bubble, and thus is irrelevant for analysing the dynamic behaviour of the interface. To perform a spectral analysis on the normal deformation
$\eta (\theta , \phi , t)$
, it is customary to expand it in a series of real spherical harmonic functions
$Y_{l}^{m}(\theta ,\phi )$
of degree
$l$
and order
$m$
, with
$0 \leqslant m \leqslant l$
:

where
$a_l^{m}(t)$
represents the time-dependent amplitude of each spherical harmonic
$Y_{l}^{m}(\theta ,\phi ) = N_{l}^{m} P_{l}^{m} (\cos (\theta )) \cos (m \phi )$
. Here
$P_{l}^{m} (\cos (\theta ))$
is a Legendre function and the normalisation factor
$N_{l}^{m} = (\max (P_{l}^{m} (\cos (\theta )) \cos (m \phi )))^{-1}$
ensures that the spherical harmonic is constrained to a maximum magnitude of unity. For a free bubble, which constitutes a closed surface, the following periodic boundary conditions must hold:


This implies that the surface deformation can only be represented by Legendre functions with integer degrees
$l$
and orders
$m$
, commonly known as associated Legendre polynomials. From the linear framework used to derive (1.1), it follows that the driving frequency and bubble radius determine the degree
$l$
of the shape mode, but not its order
$m$
. Consequently, for a given
$l$
, there are
$l+1$
possible values of
$m$
, each occurring with equal probability (see figure 1). Modes that share the same resonance frequency, and thus the same
$l$
, are referred to as degenerate. Nonlinear interactions dictate the combination of degenerate modes that ultimately emerges for each
$l$
. The first seven modes with integer degree
$l$
for
$m=0$
are depicted in figure 20(a,b). The mode spectrum of a free bubble for the first seven degrees is depicted in figure 5(a). The degree
$l$
can be interpreted as the ‘energy level’ of the mode, while the order
$m$
can be seen as the number of ‘states’ that become available with increasing
$l$
. The ‘energy level’
$l$
can be augmented by increasing either the driving frequency or the bubble radius.

Figure 5. Spectrum of allowed shape modes for the first seven
$k$
values, from 0 to 6, for (a) a free bubble, (b) a pinned wall-attached bubble, (c) a fixed-contact-angle wall-attached bubble and (d) an unconstrained wall-attached bubble. The static contact angle considered for the wall-attached bubbles is
$\alpha _0 = 3\pi /4$
. Parameters
$l$
and
$m$
are the degree and order, respectively, of the spherical harmonic
$Y_{l_{k}}^{m}$
, and
$k$
is an index that sequentially orders shape modes of a specific order
$m$
. For a free bubble, both
$l$
and
$m$
are integers, resulting in a degenerate spectrum where different shape modes share the same
$l$
. In contrast, for a pinned or fixed-contact-angle wall-attached bubble,
$l$
is generally a non-integer, leading to a non-degenerate spectrum. For an unconstrained wall-attached bubble, the spectrum is no longer quantised in
$l$
but is continuous and remains degenerate.
In general, spherical harmonics with an integer degree are not suitable for describing the surface deformation of a bubble in contact with a substrate, as they may not comply with the variational wetting condition (3.4). For a pinned bubble, the deformation basis functions must satisfy the following boundary condition:

This condition is fulfilled by Legendre functions
$P_{l}^{m} (\cos (\theta ))$
of real, though not necessarily integer, degree
$l$
that are zero at
$\theta =\alpha _0$
. Since the distinct real values of
$l$
for which this condition holds depend separately on
$m$
, they will be denoted by
$l_{k}(m)$
, where
$k$
is an integer index that sequentially orders the various roots
$l$
at a given
$m$
. The index
$k$
is chosen to start at
$m$
, analogous to the case of integer
$l$
. For a free bubble,
$k$
corresponds to
$l$
. In the azimuthal direction, periodic boundary conditions must hold:

which restrict the value of
$m$
to an integer. The first seven modes with real degrees
$l$
for
$m=0$
that comply with the condition (3.8) for
$\alpha _0 = 3\pi /4$
are depicted in figure 20(c,d). The contact angle
$\alpha _0 = 3\pi /4$
is chosen to align with the experimental observations discussed in the following sections. The mode spectrum of a pinned bubble with
$\alpha _0 = 3\pi /4$
for the first seven
$k$
values is reported in table 1(a) and also depicted in figure 5(b). In this case, shape modes are non-degenerate as each mode corresponds to a unique
$l$
value. This leads to a richer, more granular spectrum where all permissible modes can distinctly manifest. It is important to recognise that a superposition of modes sharing the same degree
$l$
but differing in order
$m$
cannot satisfy the boundary condition along the entire contact line, since
$\phi$
-dependent deformation components cannot cancel out everywhere. Only individual modes that produce uniform deformation along the contact line can meet this constraint. Compared with a free bubble, the
$l$
value for a given
$k$
value is higher, as the same number of shape-mode features must be compressed into a shorter arc of the circumference. Consequently, for a given radius of curvature, the value of
$l$
increases as the contact angle decreases. Finally, it can be observed that the deviation of
$l$
from
$k$
increases as
$k$
increases and
$m$
decreases.
If the contact line is free to move but the contact angle is fixed (
$\delta \alpha =0$
), the deformation basis functions must satisfy the following boundary condition:

Again, this condition is met by Legendre functions
$P_{l}^{m} (\cos (\theta ))$
of real, though not necessarily integer, degree
$l$
and integer order
$m$
. The first seven modes with real degrees
$l$
for
$m=0$
that comply with the condition (3.10) for
$\alpha _0 = 3\pi /4$
are depicted in figure 20(e, f). The mode spectrum of a fixed-contact-angle bubble with
$\alpha _0 = 3\pi /4$
for the first seven
$k$
values is reported in table 1(b) and also depicted in figure 5(c). The
$l$
value for a given
$k$
and
$m$
values is lower compared with the pinned-bubble case. These findings indicate that the pinning conditions can markedly influence the resonance frequency of shape modes, including the breathing mode.
In the family of Legendre functions describing the deformation of a free bubble, characterised by integer indices
$l$
and
$m$
, orthogonality is satisfied for both fixed
$l$
and fixed
$m$
:


The corresponding spherical harmonics
$Y_{l}^m = N_{l}^m P_{l}^m (\cos (\theta )) \cos (m \phi )$
are thus orthogonal to each other:

Orthogonality is also satisfied in the families of Legendre functions describing the deformation of a pinned bubble and a fixed-contact-angle bubble. However, since no function shares the same
$l$
, orthogonality holds only for fixed
$m$
:

In these families as well, the corresponding spherical harmonics are therefore orthogonal to each other:

Spherical harmonics belonging to different families are in general not orthogonal, and the inner product between them is non-zero. Consequently, the spherical harmonics in these two families are not volume-preserving as they are not orthogonal to
$Y_0^0=1$
, and therefore their 1-norm is non-zero.
If both the contact line and contact angle are free to move, the deformation can adopt the shape of a Legendre function with any real
$l$
value. This is because there are no boundary conditions imposed along the
$\theta$
direction, including the periodicity condition required for a free bubble. Consequently, as shown in figure 5(d), the mode spectrum for a bubble with a free contact line and contact angle is no longer quantised in
$l$
, but, rather, it becomes continuous. However, it remains quantised in
$m$
due to the
$2\pi$
-periodicity in the
$\phi$
direction that still holds. In this case, the spectrum exhibits degeneracy, as shape modes can share the same degree
$l$
, a scenario analogous to that of a free bubble.
4. Experimental observation of the bubble response
We present a representative example of the response of an air bubble resting against the glass substrate subjected to ultrasound driving in figure 6 (and supplementary movie 1 available at [URL]). In the figure, the top rows offer a top-view perspective captured with visible light, while the bottom rows show a side view obtained using X-ray synchrotron light. The frames have been post-processed to enhance visual clarity: noise has been reduced, and the static background has been subtracted. In the top-view images, slight shadowing effects are visible due to the angled illumination employed during capture. In the side-view images, a faint horizontal line appears as an X-ray imaging artefact, most likely resulting from the substrate not being perfectly parallel to the X-ray beam. This artefact is not related to the physical structure or properties of the bubble and does not influence the interpretation of the data presented. The equilibrium radius of the bubble is
$R_0 = {82}\;\unicode{x03BC} {\textrm{m}}$
and the equilibrium contact angle is
$\alpha _0 \approx 140^\circ$
, as directly estimated from the images. This measured contact angle is consistent with values reported in previous experimental studies involving air bubbles in water in contact with untreated glass surfaces (Ozkan & Erbil Reference Ozkan and Erbil2017). The ultrasound frequency is
$f_{ { d}} = {30}\,\textrm {kHz}$
, its peak pressure is
$p_{ { a}} = {2.7}\,{\textrm{kPa}}$
and its direction in the images is from bottom to top. The wavelength
$\lambda$
of the acoustic forcing is considerably larger than the bubble size:
$\lambda /2R_0 \sim O(10^2)$
. Therefore, the pressure field can be considered uniform around the bubble. Here, and in all subsequent figures, the time
$t = 0$
marks the onset of ultrasound excitation.

Figure 6. Overview of the response of an air bubble in contact with a glass substrate and subjected to ultrasound. The equilibrium radius of the bubble is
$R_0 = {82}\;\unicode{x03BC} {\textrm{m}}$
and the ultrasound frequency is
$f_{ { d}} = {30}\,\textrm {kHz}$
, while its peak pressure is
$p_{ { a}} = {2.7}\,{\textrm{kPa}}$
. The response can be categorised into four sequential regimes: (i) breathing mode, (ii) harmonic zonal deformation mode, (iii) half-harmonic zonal deformation mode and (iv) combined half-harmonic zonal and sectoral deformation mode. The top rows display the bubble dynamics from a top-view perspective captured with visible light, while the bottom rows present a side-view perspective, captured using X-ray synchrotron radiation. When the amplitude of the deformation modes exceeds a certain threshold, jetting can occur. The images have been denoised and their backgrounds have been removed for visual clarity.
The time evolution of the bubble response can be categorised into four distinct phases, as elaborated below and in the following section. In the initial acoustic cycles, the bubble response exclusively exhibits a breathing mode, as illustrated in the first segment of the figure. During this phase, the bubble undergoes compression and expansion cycles at the frequency of the ultrasound driving while maintaining its spherical shape. As time proceeds, the bubble tends to flatten slightly against the substrate due to the primary radiation force exerted by the ultrasound. Moreover, if the breathing mode exceeds a certain threshold in amplitude, the bubble interface destabilises leading to the development of a standing-wave pattern. Initially, this pattern is axisymmetric, its amplitude small and it matches the frequency of the acoustic driving, as outlined in the second segment of the figure. Therefore, following the nomenclature for spherical harmonics, this initial pattern is zonal and its frequency classifies it as harmonic. Then, the deformation pattern lowers its frequency to half the acoustic driving frequency and its amplitude increases, which is depicted in the third segment of the figure. The zonal standing-wave pattern is thus now half-harmonic, which is a signature trait of the Faraday instability. Interestingly, owing to the more intense shape deformation following the frequency shift, the bubble tip undergoes a vigorous folding when the shape mode reverses, ultimately resulting in the formation of a microjet that impacts the substrate, as highlighted by the red arrow in the corresponding frame. The occurrence of jetting is cyclical and does not result in bubble disruption. As the microjet manifestation is linked to the shape-mode cyclicity, its frequency is half-harmonic. The jet retraction may result in the pinching off of one or more droplets due to the Rayleigh–Plateau instability. Over time, the bubble begins to show azimuthal perturbations that break its axial symmetry and lead to the emergence of a sectoral standing-wave pattern. The zonal pattern persists but diminishes in intensity, leading to the cessation of jetting as the tip folding weakens compared with before. The mixing of a zonal pattern and a sectoral pattern marks the ultimate phase of bubble response, as reported in the fourth segment of the figure. Before the jetting stops, the bubble under examination exhibits six consecutive jets. Finally, the images further reveal that, for the glass substrate being studied, both the movement of the contact line and changes in the contact angle are allowed as the bubble interface evolves. This behaviour aligns more closely with the fully-free dynamical wetting conditions occurring in the ideal case of an unconstrained wall-attached bubble, as depicted in figure 5(d).
5. Detailed analysis of the bubble response
5.1. Breathing mode
The bubble deformation
$\eta (\theta ,\phi , t)$
can be decomposed into a volume-variant component
$\eta _{ { volume}}(\theta ,\phi , t)$
and a volume-invariant component
$\eta _{ { shape}}(\theta ,\phi , t)$
:

For a free bubble,
$\eta _{{ volume}}(\theta ,\phi , t)=a_0^0(t)Y_0^0(\theta ,\phi )$
, which corresponds to the breathing mode. The shape deformation, on the other hand, is described by
$\eta _{{ shape}}(\theta ,\phi , t)=\sum a_l^m(t) Y_l^m(\theta ,\phi )$
with
$l\neq 0$
. For a wall-attached bubble, performing such a decomposition in terms of spherical harmonics is not straightforward. As mentioned previously, there are no spherical harmonics that simultaneously satisfy the wetting condition and are volume-invariant. Nevertheless, we can still experimentally track the temporal evolution of the volumetric change, define an equivalent breathing mode amplitude and compare it with classical theoretical models for spherical bubbles. This proves useful to elucidate the accuracy of these models in the context of wall-attached bubbles and serves as a valuable means to validate the value of the acoustic driving pressure measured by the hydrophone. For this purpose, we define an average deformation of the bubble surface, analogous to an increase
$\delta R(t)$
in the radius of the spherical cap,
$ \delta R(t) = ({1}/{(2\pi (1-\cos (\alpha _0)))})\int _0^{2\pi } \int _0^{\alpha _0} \eta (\theta ,\phi ,t)\sin {\theta }\,{\textrm d}\theta \,{\textrm d}\phi .$
The volume of this equivalent spherical cap is consequently given by
$ V(t) = (({2}/{3}) - \cos (\alpha _0) + ({1}/{3})\cos ^3(\alpha _0) )\pi (R_0 + \delta R(t) )^3.$
The radius
$R_{\textit{eq}}(t)$
of a spherical bubble, which has the same volume as the equivalent spherical cap, is given by
$ R_{\textit{eq}}(t) = (({3V(t)}/{4 \pi }) ) ^{{1}/{3}} = ( ({1}/{2}) - ({3}/{4})\cos (\alpha _0) + ({1}/{4})\cos ^3(\alpha _0) )^{{1}/{3}} (R_0+ \delta R(t) ).$
For a contact angle
$\alpha _0 \approx 140^\circ$
,
$R_{\textit{eq}}(t) \approx 0.99 (R_0+ \delta R(t) )$
. Therefore, from now on, the difference between the two will be regarded as negligible, and the term
$R(t)$
will be used.
To experimentally measure
$R(t)$
we employ an image-processing algorithm to extract the bubble contour from the side-view video recordings. As long as the bubble maintains axisymmetry, i.e. before developing sectoral shape modes, its volume
$V(t)$
can be determined by integrating the area of disks stacked along the axis of symmetry
$z$
, defined by the bubble contour
$R(z,t)$
, as
$ V(t) =\pi \int _{z^{-}}^{z^{+}} R(z,t)^2 {\textrm d}z,$
where
$z^{-}$
and
$z^{+}$
represent the limiting values within which the contour of the bubble
$R(z,t)$
is defined. The radius
$R(t)$
of a volume-equivalent spherical bubble can then be determined as

The radial motion of this equivalent spherical bubble of radius
$R(t)$
in proximity to a solid boundary and subjected to a spatially uniform acoustic field can be approximated by the Rayleigh–Plesset equation for mildly compressible Newtonian fluids. The influence of the nearby solid boundary is incorporated using the method of images (Strasberg Reference Strasberg1953), yielding the following formulation:

where overdots denote time differentiation,
$\rho _{ { l}}$
the liquid density,
$c_{ { l}}$
the speed of sound in the medium,
$p_{ { g}}$
the gas pressure within the bubble,
$\sigma _{{ lg}}$
the liquid–gas surface tension,
$p_{\infty }$
the undisturbed ambient pressure,
$p_{{ d}}(t)$
the acoustic driving pressure,
${\,{\,\rm \mu}} _{ { l}}$
the dynamic viscosity of the medium and
$d=R_0$
the distance from the centre of the bubble to the wall. A full derivation of the pressure term accounting for the rigid boundary is provided in Appendix D. The method of images is typically valid when the distance between the bubble and the wall is much larger than the bubble radius. In our case, where the bubble is attached to the wall, finite-size and multipole effects may become significant. More advanced models have been developed to describe bubble dynamics near boundaries (Doinikov, Aired & Bouakaz Reference Doinikov, Aired and Bouakaz2011; Hay et al. Reference Hay, Ilinskii, Zabolotskaya and Hamilton2013), which capture more complex interactions than the simple image method. Nevertheless, despite its simplicity, the method of images accurately predicts the experimentally observed resonance size and damping of wall-attached bubbles, showing strong agreement with both radial dynamics (figure 7) and shape-mode pressure inception (figure 11), which we examine in detail later. To solve equation (5.3), it is necessary to determine the gas pressure
$p_{ { g}}$
within the bubble. Given its simplicity, the assumption of a polytropic process is widely employed:
$p_{ { g}} = p_{ { g,0}} ({R_0}/{R} )^{3n}, \ p_{ { g,0}} = p_{\infty } + ({2\sigma _{{ lg}}}/{R_0}),$
with
$p_{ { g,0}}$
denoting the bubble internal pressure at rest,
$R_0$
the equilibrium bubble radius and
$n$
the polytropic index.

Figure 7. Experimental and predicted time evolution of the equivalent radius of a bubble in contact with a glass substrate subjected to an ultrasound driving at a frequency
$f_{ { d}} = {30}\,\textrm {kHz}$
. In (a), the bubble has an an equilibrium radius
$R_0 = {89.0}\;\unicode{x03BC} {\textrm{m}}$
and the ultrasound pressure amplitude is
$p_{ { a}} = {2.1}\,{\textrm{kPa}}$
. In (b), the bubble has an an equilibrium radius
$R_0 = {102.6}\;\unicode{x03BC} {\textrm{m}}$
and the ultrasound pressure amplitude is
$p_{ { a}} = {4.6}\,{\textrm{kPa}}$
. Theoretical predictions are obtained by solving the Rayleigh–Plesset equation using the following physical parameters:
$p_{\infty } \,{=}\,{101.0}\,{\textrm{kPa}}$
,
$\sigma _{{ lg}} \,{=}\, {72}\,\textrm {mN}\,{\textrm {m}^{-1}}$
,
$c_{{ l}} \,{=}\, {1481}\,\textrm {m}\,{\textrm {s}^{-1}}$
,
${\,{\,\rm \mu}} _{ { l}} \,{=}\, {9.54}\times {10^{-4}}\,{\textrm{Pa}}\,{\textrm{s}}$
and
$\rho _{ { l}} \,{=}\, {997.8}\,\textrm {kg}\,{\textrm {m}^{-3}}$
. The prediction obtained using the polytropic gas model, with
$n = \gamma = 1.4$
, is shown by the dotted grey line. The prediction obtained using the Zhou gas model, with parameters
$\gamma =1.4$
,
$K_{{ g}} = {0.026}\,\textrm {W}\,{\textrm {m}^{-1}}{\textrm {K}^{-1}}$
,
${\mathcal{R}} = {287}\,\textrm {J}\,{\textrm {kg}^{-1}}\,{\textrm {K}^{-1}}$
and
$T_{ { g}}|_{R} = {295}\,\textrm {K}$
, is represented by the continuous black line.
In figure 7, the predictions of the theoretical model are compared with the experimental radial time evolution of bubbles subjected to ultrasound forcing at a frequency of
$f_{ { d}} = {30}\,\textrm {kHz}$
. We assume for the polytropic gas model that the gas within the bubble undergoes an adiabatic process, characterised by an exponent
$n = \gamma$
, where
$\gamma = 1.4$
denotes the specific heat ratio for the gas. This assumption is justified by the high Péclet number,
$ Pe = {R_0^2 \omega _{ { d}}}/{D_{ { g}}} \sim 100$
, where
$\omega _{ { d}}$
is the driving angular frequency and
$D_{ { g}}$
represents the thermal diffusivity of the gas. Figure 7(a) illustrates the case of a resonant-sized bubble with an equilibrium radius
$R_0 = {89.0}\;\unicode{x03BC} {\textrm{m}}$
and ultrasound pressure amplitude
$p_{{ a}} = {2.1}\,{\textrm{kPa}}$
. Here, the theoretical prediction based on the polytropic gas model significantly overestimates the experimental results and manifests a continuous growth behaviour. Conversely, figure 7(b) shows results for an over-resonant bubble with an equilibrium radius
$R_0 = {102.6}{\;\rm \unicode{x03BC}} {\textrm{m}}$
and ultrasound pressure amplitude
$p_{ { a}} = {4.6}\,{\textrm{kPa}}$
. In this case, the theoretical curve obtained using the polytropic gas model exhibits a pronounced beating pattern, which is absent in the experimental data. A beating phenomenon arises from the interference between two signals with different frequencies,
$f_1$
and
$f_2$
, where the beat frequency is given by
$f_{ { beat}} = |f_2 - f_1|$
. The observed beat frequency of approximately
${4}\,\textrm {kHz}$
suggests that the second signal, which overlaps with the forced response of the bubble at
${30}\,\textrm {kHz}$
, corresponds to the free response at the resonance frequency, which for this bubble size is approximately
${26}\,\textrm {kHz}$
. From this, we deduce that the strong overestimation seen in figure 7(a) results from constructive interference between the forced and free responses, both occurring at approximately the same frequency for this bubble size. We attribute the anomalous, underdamped free response, absent in the experimental observations, to the lack of thermal damping in the polytropic gas model. Thermal damping, as discussed by Chapman & Plesset (Reference Chapman and Plesset1971), is the dominant damping mechanism for bubbles of the sizes considered in this study.
To address the thermal interaction, we employ the Zhou model (Zhou Reference Zhou2021) for the bubble gas pressure
$p_{ { g}}$
. A detailed description of the model is provided in Appendix E. Figure 7 illustrates that the theoretical predictions derived from this model align remarkably well with the experimental results, showing a much milder beating pattern.
5.2. Meniscus waves
Upon bubble oscillations, travelling surface waves are emitted from the pinning line. These waves move towards the bubble symmetry axis,
$z$
, and equilibrate to form an axisymmetric harmonic standing wave that oscillates at the same frequency as the ultrasound applied. Figure 8 displays the initial emission of the wave, its propagation and the formation of a steady-state pattern. This phenomenon has previously been observed and investigated on flat liquid surfaces within containers, where a meniscus forms along the walls (Douady Reference Douady1990; Torres et al. Reference Torres, Pastor, Jiménez and De Espinosa1995; Batson, Zoueshtiagh & Narayanan Reference Batson, Zoueshtiagh and Narayanan2013; Ward, Zoueshtiagh & Narayanan Reference Ward, Zoueshtiagh and Narayanan2019; Shao et al. Reference Shao, Wilson, Saylor and Bostwick2021). When the container is subjected to harmonic vertical acceleration, travelling edge waves are emitted from the walls. These waves are not generated by a parametric instability but rather by synchronous adjustments of the meniscus profile in response to the harmonic modulation of the acceleration field, resulting in the harmonic emission of an axisymmetric wave from the wall in order to conserve fluid mass. Their occurrence requires no minimum threshold in the magnitude of the driving acceleration, and their amplitude is linearly proportional to this magnitude. The meniscus forms when the contact angle deviates from
$90^\circ$
, and meniscus waves occur regardless of whether the contact line is pinned or free. These waves can be suppressed by filling the liquid container to the brim, thereby pinning the contact line at the brim with a contact angle of
$90^\circ$
. To our knowledge, our study is the first to recognise the presence of meniscus waves on ultrasound-driven bubbles in contact with a solid surface. In this case, the presence of a wall mandates that the acceleration of the bubble interface near the wall must be parallel to the wall, as dictated by the non-penetration condition. As a result, when the contact angle deviates from
$90^\circ$
, the acceleration along the pinning line of the bubble is not orthogonal to the bubble surface and meniscus waves are emitted. Conversely, when a bubble in contact with a wall exhibits a contact angle of
$90^\circ$
meniscus waves are not expected to occur, as the direction of acceleration is orthogonal to the bubble surface at all points.

Figure 8. Time evolution of harmonic meniscus waves produced by the harmonic acceleration of the bubble surface close to the wall. The wave pattern is highlighted by performing a subtraction operation between the X-ray image capturing the wave at its peak and that with the wave at its trough. (a) Initial emission of the travelling meniscus wave. (b) Propagation of the travelling meniscus wave. (c) Steady state of meniscus waves forming a standing-wave pattern.
5.3. Faraday waves
When the amplitude of the ultrasound pressure exceeds a critical threshold, a half-harmonic parametric instability known as Faraday instability emerges, supplanting the harmonic pattern of meniscus waves. The resulting standing-wave pattern, also termed shape mode, is the result of the nonlinear interaction between the external periodic forcing and the natural frequencies of the bubble interface. This interaction leads to the amplification of certain modes while damping others, resulting in the observable patterns. Initially, the Faraday wave pattern exhibits axisymmetry, akin to meniscus waves, resulting in a purely zonal shape mode. The wavelength of Faraday waves is approximately twice that of meniscus waves, and its amplitude is substantially higher, as evidenced in figure 9(a,b). As time proceeds, Faraday waves progressively lose their axisymmetry, ultimately stabilising into a non-axisymmetric shape mode. This finding supports the theoretical prediction for spherical entities as stated by Chossat et al. (Reference Chossat, Lauterbach and Melbourne1991), which asserts that no axisymmetric mode with
$l \geqslant 2$
is inherently stable. The ultimate shape mode of the bubble results from the superposition of a sectoral mode onto the pre-existing zonal mode, as depicted in figure 9(c). This transition is evidenced by the development of lobes in the azimuthal direction, which overlap with the existing axisymmetric lobes in the polar direction. The absence of layer-symmetric patterns indicates that tesseral shape modes are not present. This contrasts with the findings of Fauconnier et al. (Reference Fauconnier, Béra and Inserra2020), who reported the presence of tesseral modes. However, their study employed a different experimental set-up: the substrate was made of PMMA, the ultrasound was oriented parallel to the substrate and shape modes were identified solely from a top-view perspective.

Figure 9. Comparison between meniscus waves and Faraday waves on the same bubble as that of figure 8. (a) Meniscus waves are harmonic waves whose emergence is not contingent upon surpassing a minimum amplitude threshold in the ultrasound driving. (b) Faraday waves are half-harmonic waves characterised by a wavelength approximately double that of meniscus waves and a significantly higher amplitude. Their occurrence is contingent upon surpassing a critical threshold in the ultrasound driving amplitude. Initially, Faraday waves are axisymmetric, i.e. the shape mode is purely zonal. (c) Over time, Faraday waves gradually lose their axisymmetry, manifesting a sectoral mode that overlays the pre-existing zonal mode. This transition is marked by the emergence of lobes in both the polar and azimuthal directions. The persistence of the bottom lobe indicates the absence of tesseral modes.
Given the truncated geometry of a wall-attached bubble, the degree
$l$
of the observed zonal shape mode, generally a real number, can be estimated only by measuring the spacing between the ridges of the shape mode as seen from the side view and comparing it with the circumference of the bubble. The length of the arc
$a$
subtended by the chord of length
$c$
joining two ridges of the shape mode is calculated as

Thus,
$l$
is given by the number of times the arc length can divide the bubble circumference:


Figure 10. Experimentally measured shape mode degree (a) and order (b) of bubbles in contact with a glass substrate, subjected to an ultrasound driving frequency of
$f_{ { d}} = {30}\,\textrm {kHz}$
, as a function of their equilibrium radius. The degree varies continuously with the bubble radius, whereas the order is quantised and only assumes integer values.
Conversely, the order
$m$
of the observed sectoral shape mode can be readily determined by counting the number of lobes along the azimuthal direction when the bubble is viewed from above. Since the bubble contour remains closed and periodic in this view,
$m$
must be an integer. At a fixed ultrasound driving frequency
$f_{ { d}}={30}\,\textrm {kHz}$
, the measured values of
$l$
vary approximately linearly with the bubble equilibrium radius
$R_0$
, while the values of
$m$
show a stepwise progression, as illustrated in figure 10. The fact that
$l$
spans a practically continuous range – unlike the quantised behaviour expected under strict boundary conditions – suggests that the contour endpoints of the bubble are effectively free. This interpretation is supported by the observed motion of the contact line and changes in contact angle during bubble oscillation. Thus, the current configuration can be considered as corresponding to the idealised scenario of a wall-attached bubble with free boundary conditions, exhibiting a spectrum of shape modes that is continuous in
$l$
and discrete in
$m$
, as discussed in § 3.2 and shown in figure 5(d). It can be determined from figure 10 that, for a given radius
$R_0$
, the value of
$m$
of the sectoral mode that manifests is approximately the integer closest to the observed value of
$l$
for that radius. Therefore, the bubble deformation
$\eta (\theta ,\phi ,t)$
can be decomposed as

where
$\lfloor \rceil$
indicates the rounding operator. It is important to note that, for each value of
$l$
, only a specific pair of values of
$m$
manifest from the set of possible values, specifically the zeroth value and the highest possible value of
$m$
for the given
$l$
. Consequently, under the present wetting conditions, the spectrum of shape modes is degenerate, similar to that of a free bubble (Cattaneo et al. Reference Cattaneo, Guerriero, Shakya, Krattiger, G. Paganella, Narciso and Supponen2025), with the difference that the values of
$l$
are real and continuous rather than integer.

Figure 11. Shape-mode order
$m$
of bubbles in contact with a glass substrate, subjected to an ultrasound driving frequency
$f_{ { d}} = {30}\,\textrm {kHz}$
, as a function of the equilibrium radius and the ultrasound pressure amplitude. Experimental data are depicted with symbols, while the theoretical threshold for the onset of shape modes is based on the model from Francescutto & Nabergoj (Reference Francescutto and Nabergoj1978), as detailed in the text. The solid lines correspond to a surface tension value
$\sigma _{{ lg}} = {65}\,\textrm {mN}\,{\textrm {m}^{-1}}$
, and the dotted lines represent
$\sigma _{{ lg}} = {72}\,\textrm {mN}\,{\textrm {m}^{-1}}$
. The other parameters used in the model are:
$\rho _{ { l}} = {997.8}\,\textrm {kg}\,{\textrm {m}^{-3}}$
,
${\,\rm \mu} _{{ l}} = {9.54}\times {10^{-4}}\,{\textrm{Pa}}\,{\textrm{s}}$
and
${\,\rm \mu} _{ { e}}$
extracted from the theoretical plots in Chapman & Plesset (Reference Chapman and Plesset1971).
The occurrence of the Faraday instability, as detected experimentally in relation to bubble size and ultrasound pressure amplitude, is summarised in figure 11. The instability region appears to be arranged into distinct tongues, each corresponding to the specific order
$m$
observed. This finding is notable as the presence of distinct tongues indicates that the shape modes exhibit discrete resonance sizes for a fixed acoustic driving frequency. Given that the accessible range of shape-mode degrees
$l$
is continuous, one would expect the associated resonant bubble sizes to be similarly continuous. However, the quantisation of the azimuthal wavenumber
$m$
introduces a selectivity on the bubble size, resulting in a discrete number of resonance sizes. The experimentally observed instability thresholds agree well with those predicted by the theory from Francescutto & Nabergoj (Reference Francescutto and Nabergoj1978) and Nabergoj & Francescutto (Reference Nabergoj and Francescutto1978) developed for shape oscillations of free bubbles. Full details of our implementation of the Francescutto & Nabergoj model for wall-attached bubbles are provided in Appendix F. The minimum inception threshold for shape modes occurs at the bubble radius that corresponds to the resonant radius, which is also close to the resonance size for the shape mode with
$m=4$
, making it the most easily excited shape mode. We note that to achieve alignment with the experimental data, we found it necessary to use a surface tension value of
${65}\,\textrm {mN}\,{\textrm {m}^{-1}}$
. This is lower than the surface tension of a clean air–water interface, which is around
${72}\,\textrm {mN}\,{\textrm {m}^{-1}}$
. The reduction in surface tension could be attributed to the presence of surfactants in the liquid medium that have adsorbed onto the bubble interface and thus decreased its surface tension. This value closely aligns with that reported in Versluis et al. (Reference Versluis, Goertz, Palanchon, Heitman, Van Der Meer, Dollet, De Jong and Lohse2010), where the authors investigated the shape modes of free air bubbles in an aqueous medium driven by ultrasound. The theory predicted by Francescutto & Nabergoj (Reference Francescutto and Nabergoj1978) achieves strong agreement with experiments even without accounting for meniscus waves, suggesting that these waves are unlikely to be a significant factor in triggering shape modes. This observation aligns with the findings of Batson et al. (Reference Batson, Zoueshtiagh and Narayanan2013), who determined in the context of flat water–air interfaces that significant influence of meniscus waves is more likely to occur for harmonic axisymmetric Faraday waves, as they share the same frequency and pattern with meniscus waves.

Figure 12. Time evolution of the experimental amplitude of the breathing, zonal and sectoral modes for the case
$l \approx m=3$
(
$l =3.22$
). The red dotted line marks the threshold amplitude of the breathing mode for the onset of shape modes. The amplitude of the breathing mode is shown as dashed during the occurrence of the sectoral mode, indicating that accurate measurement is not possible due to the loss of bubble axisymmetry. Below is a sequence of top-down optical and side-view X-ray images illustrating the dynamics of the bubble. The images have been denoised and their backgrounds have been removed for visual clarity. The equilibrium bubble radius is
$R_0 = {72} {\;\rm \unicode{x03BC}} {\textrm{m}}$
. The ultrasound driving frequency is
$f_{{ d}} = {30}\,\textrm {kHz}$
and its amplitude is
$p_{ { a}} = {5.3} \,{\textrm{kPa}}$
.

Figure 13. Time evolution of the experimental amplitude of the breathing, zonal and sectoral modes for the case
$l \approx m=4$
(
$l = 4.21$
). The red dotted line marks the threshold amplitude of the breathing mode for the onset of shape modes. The amplitude of the breathing mode is shown as dashed during the occurrence of the sectoral mode, indicating that accurate measurement is not possible due to the loss of bubble axisymmetry. Below is a sequence of top-down optical and side-view X-ray images illustrating the dynamics of the bubble. The images have been denoised and their backgrounds have been removed for visual clarity. The equilibrium bubble radius is
$R_0 = {89} {\;\rm \unicode{x03BC}} {\textrm{m}}$
. The ultrasound driving frequency is
$f_{ { d}} = {30}\,\textrm {kHz}$
and its amplitude is
$p_{ { a}} = {2.1} \,{\textrm{kPa}}$
.

Figure 14. Time evolution of the experimental amplitude of the breathing, zonal and sectoral modes for the case
$l \approx m=5$
(
$l =4.85$
). The red dotted line marks the threshold amplitude of the breathing mode for the onset of shape modes. The amplitude of the breathing mode is shown as dashed during the occurrence of the sectoral mode, indicating that accurate measurement is not possible due to the loss of bubble axisymmetry. Below is a sequence of top-down optical and side-view X-ray images illustrating the dynamics of the bubble. The images have been denoised and their backgrounds have been removed for visual clarity. The equilibrium bubble radius is
$R_0 = {103} {\;\rm \unicode{x03BC}} {\textrm{m}}$
. The ultrasound driving frequency is
$f_{{ d}} = {30}\,\textrm {kHz}$
and its amplitude is
$p_{ { a}} = {4.4} \,{\textrm{kPa}}$
.
5.4. Time evolution of shape modes
The time evolution of the amplitude of bubble deformation components, specifically the breathing, zonal and sectoral modes, is presented in figures 12–15 for cases with
$l \approx m = 3$
,
$l \approx m = 4$
,
$l \approx m = 5$
and
$l \approx m = 6$
, respectively. The amplitude of the breathing mode
$a_0^0(t)$
is approximately measured using the method outlined in § 5.1. This approach remains accurate as long as the bubble maintains an approximately axisymmetric shape. The amplitude
$a_l^0(t)$
of a zonal shape mode
$Y_{l}^{0}(\theta ,\phi )=P_{l}^{0}(\cos (\theta ))$
can be determined by calculating the inner product between the cross-section of the bubble deformation,
$\eta (\theta ,0,t)$
, along the
$z$
axis, and the shape mode itself:

where
$ {(2l+1)}/{2}$
is a normalisation factor that follows from Rodrigues’ formula. The amplitude
$a_m^m(t)$
of a sectoral shape mode
$Y_{m}^{m}(\theta ,\phi )= N_{m}^{m} P_{m}^{m} (\cos (\theta )) \cos (m \phi )$
can be computed in a similar way by computing the inner product between the cross-sections along the
$x$
axis and
$y$
axis of the bubble deformation
$\eta (\pi /2,\phi ,t)$
and of the shape mode
$Y_{m}^{m}(\pi /2,\phi )=\cos (m\phi )$
, respectively, as follows:

However, under the present conditions, quantifying the amplitude of the zonal mode using (5.7) is unfeasible. The non-quantised nature of the shape mode precludes precise measurement of its degree. Moreover, the bubble deformation, constrained within the range
$\theta = \pm \alpha _0$
, renders the calculation inherently ill-posed. Finally, the concurrent presence of a sectoral shape mode complicates the isolation of its specific contribution to the cross-sectional deformation. To address these challenges, we estimate the amplitude of the shape mode by measuring the amplitude of the lower-lobe oscillation, at
$\theta = 0$
. By definition, it directly corresponds to the amplitude of the shape mode, and it remains unaffected by sectoral modes, as these modes exhibit zero deformation at
$\theta = 0$
. These issues do not arise in estimating the sectoral mode amplitude due to the quantised nature of the shape-mode order, the periodicity of the deformation and the fact that the zonal mode merely introduces a uniform offset to the deformation in the examined cross-section.

Figure 15. Time evolution of the experimental amplitude of the breathing, zonal and sectoral modes for the case
$l \approx m=6$
(
$l =6.02$
). The red dotted line marks the threshold amplitude of the breathing mode for the onset of shape modes. The amplitude of the breathing mode is shown as dashed during the occurrence of the sectoral mode, indicating that accurate measurement is not possible due to the loss of bubble axisymmetry. Below is a sequence of top-down optical and side-view X-ray images illustrating the dynamics of the bubble. The images have been denoised and their backgrounds have been removed for visual clarity. The equilibrium bubble radius is
$R_0 = {129}{\;\rm \unicode{x03BC}} {\textrm{m}}$
. The ultrasound driving frequency is
$f_{{ d}} = {30}\,\textrm {kHz}$
and its amplitude is
$p_{ { a}} = {10.5} \,{\textrm{kPa}}$
.
In all examined cases, the amplitude of the breathing mode ranges between 1.5 and 2.4 times the threshold required to trigger shape modes for the specific bubble size under investigation. The temporal evolution of shape modes exhibits consistent patterns across all cases. The zonal mode invariably emerges first, followed by the sectoral mode. At the breathing mode amplitudes relative to the shape-mode threshold employed, the zonal mode typically occurs after 10–14 ultrasound cycles. A notable exception is observed in the case of
$l \approx m = 3$
, where the zonal mode initiates after just 4 cycles. Following the onset of the zonal mode, the sectoral mode typically manifests after 6–10 ultrasound cycles. The amplitude of the zonal mode increases approximately linearly with time until the inception of the sectoral mode. At its peak, the zonal mode is 6–8 times more intense than the breathing mode. This factor appears to increase as the interval between the onset of the zonal mode and the initiation of the sectoral mode prolongs. The pronounced amplification of bubble deformation produced by shape modes arises from their intrinsic capacity to concentrate kinetic energy at specific regions on the bubble surface, the shape-mode lobes. This contrasts with the breathing mode, where kinetic energy and deformation are distributed uniformly across the entire bubble surface. The concentration of energy at specific points produced by shape modes is the key factor that enables jet formation at low-amplitude ultrasound pressures, as shown in this study. Upon the onset of the sectoral mode, the amplitude of the zonal mode ceases to increase and begins to decline, stabilising at approximately 0.5–0.7 times its prior peak value. Simultaneously, the intensity of the sectoral mode grows to a level comparable to that of the zonal mode, indicating a possible energy transfer from the zonal to the sectoral mode. This decrease in the overall shape-mode intensity suggests that the most favourable conditions for jet formation occur when only the zonal mode is active, as this configuration results in the highest deformation amplitude. For the remainder of the recording period, the zonal and sectoral modes exhibit minimal variation, maintaining a largely stable pattern. We do not observe significant changes in the amplitude of the breathing mode when the zonal mode develops. Once the sectoral mode emerges, extracting volume becomes troublesome, as the bubble shape loses axisymmetry. Therefore, we can only assume that the amplitude remains constant during this regime as well.
The delayed emergence of the ultimate shape-mode pattern observed in wall-bounded bubbles contrasts with the behaviour of free bubbles, where the final Faraday pattern appears directly, without transitioning through an intermediate stage of purely axisymmetric shape modes (Cattaneo et al. Reference Cattaneo, Guerriero, Shakya, Krattiger, G. Paganella, Narciso and Supponen2025). One possible explanation is the friction at the contact line, which may restrict the bubble oscillatory motion along the azimuthal direction, thus hindering the onset of surface destabilisation in that direction. In contrast, motion along the polar direction is not similarly constrained, making zonal modes easier to develop initially. These zonal modes, superimposed on the breathing-mode spherical oscillation, may then generate sufficient surface deformation to eventually trigger the onset of sectoral modes. Interestingly, this progression is not unprecedented – it has been observed in sessile droplets on vibrating surfaces, where zonal modes appear first. Subsequently, the radial oscillation of these axisymmetric waves induces azimuthal waves at the contact line, leading to the formation of sectoral patterns (Vukasinovic et al. Reference Vukasinovic, Smith and Glezer2007; Panda et al. Reference Panda, Kahouadji, Tuckerman, Shin, Chergui, Juric and Matar2023).
5.5. Instability threshold for repeated jets
When the amplitude of the shape mode exceeds a critical threshold, the inward folding of the bottom shape-mode lobe generates a high-velocity jet directed towards the substrate. After onset, jetting occurs at each cycle of the shape mode. Given that the shape mode exhibits a half-harmonic behaviour, the jet frequency is half that of the applied ultrasound driving. Shape-mode-induced cyclic jets on bubbles have likely been observed in previous studies that documented repeated jets resulting from surface deformations at driving frequencies spanning hertz (Crum Reference Crum1979), kilohertz (Prabowo & Ohl Reference Prabowo and Ohl2011) and megahertz (Vos et al. Reference Vos, Dollet, Versluis and De Jong2011) ranges.

Figure 16. (a) X-ray image sequence of shape-mode-induced jetting. The bottom lobe undergoes maximum excursion, and when it rapidly folds inwards, a singularity forms at the fluid interface, where both surface curvature and velocity diverge, leading to the formation of a high-speed jet directed towards the substrate, driven by the self-focusing of kinetic energy at the surface singularity. The equilibrium radius of the bubble is
${121}{\;\rm \unicode{x03BC}} {\textrm{m}}$
, the driving pressure is
$p_{ { a}} = {9.4}\,{\textrm{kPa}}$
and the frequency is
$f_{ { d}} = {30}\,\textrm {kHz}$
. (b) X-ray image sequence of ultrasound-driven inertial jetting. The bubble forms a jet as soon as the ultrasound amplitude reaches its steady state, leaving no time for shape modes to develop. This produces a single, quick jet followed by rapid bubble fragmentation. The equilibrium radius of the bubble is
${91}{\;\rm \unicode{x03BC}} {\textrm{m}}$
, the driving pressure is
$p_{ { a}} = {17}\,{\textrm{kPa}}$
and the frequency is
$f_{ { d}} = {30}\,\textrm {kHz}$
. The background of the images has been removed for visual clarity.
Figure 16(a) illustrates that, for a bubble with an equilibrium radius
$R_0 = {121}{\;\rm \unicode{x03BC}} {\textrm{m}}$
, subjected to a driving pressure
$p_{ { a}} = {9.4}\,{\textrm{kPa}}$
and frequency
$f_{ { d}} = {30}\,\textrm {kHz}$
, just before jet emission occurs, the bottom shape-mode lobe folds inward, deforming into an approximate truncated cone. This creates a singularity point at the fluid interface, characterised by divergent surface curvature and velocity. The resulting self-focusing of kinetic energy at this singularity point initiates the jet formation. Jets originating from the bottom shape-mode lobe are predominantly observed due to its higher amplitude compared with other lobes, as it is not constrained by the substrate. This contrasts with free bubbles, where all lobes exhibit equal amplitude, allowing simultaneous jet formation when their amplitude is sufficiently high (Cattaneo et al. Reference Cattaneo, Guerriero, Shakya, Krattiger, G. Paganella, Narciso and Supponen2025). This jetting phenomenon can be considered the spherical analogue of the jets formed by Faraday waves on vertically vibrating liquid baths (Longuet-Higgins Reference Longuet-Higgins1983; Zeff et al. Reference Zeff, Kleber, Fineberg and Lathrop2000). Jets arising from a collapsing conical interface are also observed in other phenomena, including bubble bursting at fluid interfaces (Kientzler & Arons Reference Kientzler and Arons1954; Duchemin et al. Reference Duchemin, Popinet, Josserand and Zaleski2002), droplet impact on liquid pools (Worthington & Cole Reference Worthington and Cole1897; Thoroddsen et al. Reference Thoroddsen, Takehara, Nguyen and Etoh2018), cavitation bubble collapse in extremely close proximity to solid boundaries (Lechner et al. Reference Lechner, Lauterborn, Koch and Mettin2019; Reuter & Ohl Reference Reuter and Ohl2021), cavitation-bubble-pair interactions (Mishra et al. Reference Mishra, Bourquard, Roy, Lakkaraju, Ghosh and Supponen2022; Fan et al. Reference Fan, Bussmann, Reuter, Bao, Adami, Gordillo, Adams and Ohl2024) and bubble coalescence (Jiang et al. Reference Jiang, Rotily, Villermaux and Wang2024).
Shape-mode-induced jets are fundamentally different from classical inertial jets, which require considerably higher acoustic pressures to form – almost an order of magnitude more. Figure 16(b) illustrates an example of ultrasound-driven inertial jetting for a bubble with an equilibrium radius
$R_0 = {91}{\;\rm \unicode{x03BC}} {\textrm{m}}$
, which is near its resonant size. The bubble is subjected to a driving pressure
$p_{ { a}} = {17}\,{\textrm{kPa}}$
and frequency
$f_{ { d}} = {30}\,\textrm {kHz}$
. Here, the bubble jets immediately upon reaching the steady-state ultrasound amplitude, leaving no time for shape modes to develop. The result is a single, swift jet followed by rapid fragmentation of the bubble. Unlike jets formed by shape modes, this jet arises from the intense pressure gradient applied to the bubble rather than from an instability at the interface. Thus, shape-mode-induced jets represent a unique mechanism for concentrating energy through relatively low acoustic excitation, marking a distinct contrast to the behaviour of classical inertial jets.

Figure 17. (a) Jetting occurrence as a function of bubble size, shape-mode order and driving ultrasound pressure. A minimum in driving pressure occurs at the resonant radius. (b) Jetting occurrence as a function of bubble size, shape-mode order and acceleration of the bottom shape-mode lobe. A distinct threshold in acceleration, regardless of bubble size, is identified.
In figure 17(a), we illustrate the relationship between bubble size and the ultrasound pressure needed to cause shape-mode-induced jet formation. The required pressure reaches a minimum at the resonant radius for the applied driving frequency. In figure 17(b), we identify a distinct threshold in the acceleration of the bottom shape-mode lobe of approximately
$0.2$
μm μs−2, irrespective of the bubble size, beyond which jetting occurs. The lobe acceleration
$\ddot z_{ { lobe}}$
is estimated from the zonal shape mode amplitude as
$\ddot z_{{ lobe}} \approx \omega _{0,l}^2 a_{l}^0$
, where
$\omega _{0,l} = {\omega _{\textrm { {d}}}}/{2}$
is the angular frequency of the shape mode. This consistent acceleration threshold arises because Faraday instability-driven jetting depends solely on the motion of individual shape-mode lobes, which is independent of bubble size. While increasing the bubble size increases the number of lobes, their individual width – set by the driving frequency – remains unchanged. As a result, the critical lobe height and corresponding acceleration required for jetting remain constant across bubble sizes. This behaviour is fundamentally distinct from that of inertial jets, where jet strength scales with the pressure difference across the bubble, and therefore depends on both the magnitude of the pressure gradient and the bubble size (Supponen et al. Reference Supponen, Obreschkow, Tinguely, Kobel, Dorsaz and Farhat2016).
The speed of jets generated by shape modes is found to reach
$u_{ { jet}} = {30}{\,\textrm {m s}^{-1}}$
. Using the fluid hammer pressure formula,
$p_{ { wh}} = \rho _{ { l}} c_{ { l}} u_{{ jet}}$
, for a water flow impinging on a rigid material, the estimated impact pressure on the substrate is calculated to be up to 45 MPa. This suggests that even at relatively low ultrasound driving pressures, the repeated jet formations are sufficiently powerful to cause damage to biological and ductile materials. This mechanism may hold significant importance in ultrasound-based cleaning applications and, as demonstrated in our previous work, plays a pivotal role in sonoporation for targeted drug delivery (Cattaneo et al. Reference Cattaneo, Guerriero, Shakya, Krattiger, G. Paganella, Narciso and Supponen2025). As it displaces, the jet may eventually break into multiple droplets following the Rayleigh–Plateau instability. Jetting persists indefinitely unless the emergence of the sectoral shape mode dampens the lobe amplitude below the critical threshold, or the bubble shape becomes too chaotic as a result of the disruptive influence of jets on the shape itself.

Figure 18. Comparison of the time evolution of the bubble interface between experimental observations and numerical simulations using the BEM, shown from both top- and side-view perspectives. (a) Transition from a purely zonal shape mode to a combined zonal and sectoral mode. The bubble has an equilibrium radius
$R_0 = {120}\,\rm {\;\rm \unicode{x03BC}} {\textrm{m}}$
and is driven by ultrasound with an amplitude
$p_{{ a}} = {9.4}\,{\textrm{kPa}}$
at a frequency
$f_{ { d}} = {30}\,\rm {kHz}$
. (b) Jetting formation during an almost purely zonal shape mode. The bubble has an equilibrium radius
$R_0 = {129}\,\rm {\;\rm \unicode{x03BC}} {\textrm{m}}$
and is driven by ultrasound with an amplitude
$p_{ { a}} = {10.5}\,{\textrm{kPa}}$
at a frequency
$f_{{ d}} = {30}\,\textrm {kHz}$
.
6. Numerical simulations
We now complement our experimental results with three-dimensional BEM numerical simulations (e.g. Wang & Manmi Reference Wang and Manmi2014), modelling the evolution of the bubble interface and formation of jets. Full details of the BEM model employed are provided in Appendix G. Figure 18 presents a sequence of phase-matched frames comparing the experimentally observed and BEM-calculated time evolution of the bubble interface. In figure 18(a), we examine the transition from a purely zonal mode to a combined zonal and sectorial mode for a bubble with an equilibrium radius of
$R_0 = {129}\,\rm {\;\rm \unicode{x03BC}} {\textrm{m}}$
, subjected to ultrasound driving at an amplitude
$p_{ { a}} = {10.5}\,{\textrm{kPa}}$
and a frequency
$f_{ { d}} = {30}\,\textrm {kHz}$
. The simulation closely matches the observed data, in terms of both the onset of the transition, occurring at around
$t={500}\rm {\;\rm \unicode{x03BC}} \textrm {s}$
(15 ultrasound cycles), and the transition duration, which spans about
${400}\rm {\;\rm \unicode{x03BC}} \textrm {s}$
. Moreover, throughout the entire transition, the amplitude of the various shape modes exhibits good agreement between the simulation and experimental results. In figure 18(b), we examine the formation of a jet during an almost purely zonal shape mode – the most common manifestation, as previously observed – for a bubble with an equilibrium radius of
$R_0 = {120}\,\rm {\;\rm \unicode{x03BC}} {\textrm{m}}$
, subjected to ultrasound driving at an amplitude
$p_{ { a}} = {9.4}\,{\textrm{kPa}}$
and a frequency
$f_{ { d}} = {30}\,\textrm {kHz}$
. The simulation terminates when a drop detaches from the first ejected jet due to the Rayleigh–Plateau instability, causing the mesh to break up. Up to that point, the simulation and experimental observations of the bubble interface evolution are in good agreement. Note that in the final time instant compared, the discrepancy between the experimental and simulation results, especially clear in the jet displacement, arises because the two frames are not phase-matched as the simulation concludes 6
${\;\rm \unicode{x03BC}}$
s earlier. Although the simulation predicts the emission of the first jet five ultrasound cycles earlier than observed experimentally, this discrepancy is acceptable given the inherent variability in the number of cycles required to produce a jet in the experiments. In conclusion, these results support the use of the BEM for simulating the dynamics of wall-attached ultrasound-driven bubbles of this scale.

Figure 19. Simulated examples of shape-mode-induced jetting, modelled using the BEM, presented in both top and side views for various bubble sizes corresponding to the following regimes: (a)
$m = 3$
, (b)
$m = 4$
, (c)
$m = 5$
and (d)
$m = 6$
. Images are shown at intervals of one per ultrasound cycle, resulting in an interframe time of 33.3
${\;\rm \unicode{x03BC}}$
s, except for the last image, which captures the last moment before numerical mesh rupture occurs following jet formation. The bubbles have equilibrium radius
$R_0$
of 70, 90, 115 and 135
${\;\rm \unicode{x03BC}}$
m, and are subjected to pressure amplitudes
$p_{{ a}}$
of 18, 3.6, 14 and 24 kPa, respectively, at a frequency
$f_{ { d}}$
of 30 kHz.
Figure 19 displays examples of simulated shape-mode-induced jetting, modelled using the BEM, for various bubble sizes corresponding to the
$m = 3$
,
$m = 4$
,
$m = 5$
and
$m = 6$
regimes. The simulation stops when the jet either touches the substrate or pinches off a droplet, disrupting the mesh and thereby restricting the simulation to display only the first jet occurrence. We select ultrasound pressures to induce shape-mode-induced jet formation within a similar number of ultrasound cycles across different cases. For the representative bubble in the
$m=4$
regime, a pressure of just
$p_{ { a}} = {3.6}\,\rm {kPa}$
over eight ultrasound cycles is sufficient to initiate jetting. This is because the bubble radius,
$R_0 ={90}\rm {\;\rm \unicode{x03BC}} {\textrm{m}}$
, nearly matches the resonant radius. In the other regimes, higher ultrasound pressures are needed to achieve jet formation within a similar cycles count:
$p_{ { a}} ={14}\,{\textrm{kPa}}$
for the
$m=5$
case (
$R_0 ={115} \;\rm {\rm \unicode{x03BC}} {\textrm{m}}$
),
$p_{ { a}} ={18}\,{\textrm{kPa}}$
for the
$m=3$
case (
$R_0 ={70}\rm \;{\rm \unicode{x03BC}} {\textrm{m}}$
) and
$p_{ { a}} ={24}\,{\textrm{kPa}}$
for the
$m=6$
case (
$R_0 ={135}\;\rm {\rm \unicode{x03BC}}{\textrm{m}}$
). For these bubble sizes, jetting can occur even at lower pressures, but initiating it requires more cycles. In any case, the pressure needed to generate jets through shape modes is significantly lower than the pressure required for classical inertial jets.
7. Conclusions
In this study, we experimentally investigated the behaviour of individual micrometric air bubbles resting against a rigid substrate when subjected to ultrasound driving, focusing specifically on the development of interfacial instabilities over time. To achieve this, we developed a dual-view imaging set-up that integrates a visible-light top view with a phase-contrast X-ray side view, powered by the synchrotron radiation source at the ESRF on beamline ID19. This system enables high-speed, high-magnification simultaneous recordings from both perspectives, providing unprecedented insights into the dynamic evolution of bubble shapes under ultrasound driving.
The observed bubble dynamics can be organised into four distinct consecutive regimes: (i) the initial purely spherical oscillation, also referred to as breathing mode, (ii) the onset of harmonic axisymmetric meniscus waves, (iii) the emergence of a higher-amplitude half-harmonic axisymmetric mode, also referred to as zonal mode, triggered by the Faraday instability and (iv) the superposition of a half-harmonic sectoral mode onto the zonal mode. Notably, the stepwise transition through distinct shape-mode regimes contrasts with the behaviour of free bubbles, which exhibit their ultimate Faraday wave pattern immediately upon the onset of instability.
When a lobe of the Faraday-induced shape mode folds inward with sufficient intensity during its oscillatory motion, it generates a jet directed towards the substrate. This jet formation occurs at mild ultrasound pressures and recurs with each cycle of the shape mode. These jets, arising from interface instabilities, represent a new and distinct class of bubble-generated jets. They fundamentally differ from classical inertial jets, which form under high-pressure gradients during the bubble collapse phase and require significantly higher acoustic pressures.
Our detailed analysis of each regime reveals several important insights. To accurately model the bubble volume evolution over time, the commonly used polytropic gas approximation in bubble dynamics is insufficient. Accurate predictions are only achieved with a more advanced model that accounts for thermal damping effects. The observed harmonic axisymmetric pattern is driven by meniscus waves originating from the contact line. These waves are equivalent to those found on flat water–air interfaces in containers with a wall-contact angle different from
$90^\circ$
, when subjected to vertical oscillations. The shape-mode spectrum is derived through a kinematic analysis under various wetting conditions. For pinned contact line or fixed contact angle boundary conditions, the spectrum is quantised and non-degenerate, whereas for free boundary conditions, the spectrum becomes degenerate and exhibits a continuous range of shape-mode degrees. Our experimental results, consistent with the non-restricting boundary conditions imposed by the chosen substrate, confirm this expected continuity in shape-mode degrees and the observed degeneracy of the shape modes. This contrasts with free bubbles, which, while also exhibiting a degenerate spectrum, display only discrete shape-mode degrees constrained by their spherical periodicity.
Furthermore, we find a strong correlation between the experimentally measured ultrasound pressure thresholds for the onset of shape modes and the theoretical predictions based on the model from Francescutto & Nabergoj (Reference Francescutto and Nabergoj1978), when using a lowered surface tension value of
${65}\,\textrm {mN}\,{\textrm {m}^{-1}}$
. This reduction in surface tension might be attributed to the presence of surfactants on the bubble surface. Our analysis of the time evolution of breathing, zonal and sectoral mode amplitudes reveals that shape modes reach significantly higher amplitudes compared with the breathing mode. This suggests that shape modes are more effective at amplifying acoustic energy and concentrating it in localised regions, which explains the formation of shape-mode-driven jets under moderate acoustic pressures. We also identify a critical acceleration threshold for shape-mode lobes, independent of bubble size, beyond which jetting occurs. Unlike in free bubbles, where jets may form from any point on the surface, jetting in wall-attached bubbles consistently emerges from the side not constrained by the substrate.
Finally, we supplement our experimental results with three-dimensional simulations of bubble dynamics using the BEM which is based on the assumption of potential flow. These simulations align closely with the experimental visualisations, accurately capturing the evolution of bubble shape, the event timing and the required ultrasound pressure, in both the transition from zonal to mixed modes (zonal and sectoral) and the formation of jets following the emergence of strong zonal modes. Our results support the use of the BEM as a reliable approach for modelling the dynamics of ultrasound-driven wall-attached bubbles at this scale.
In conclusion, our study provides valuable insights into the dynamics of wall-attached bubbles under ultrasound driving, advancing our understanding of the development of surface instabilities and jet formation. This research has significant implications for industrial and biomedical applications, such as ultrasound cleaning and bubble-assisted targeted drug delivery.
Supplementary movie
The supplementary movie is available at https://doi.org/10.1017/jfm.2025.10457.
Acknowledgements
The results presented in this work were obtained during the allocated beamtime for proposal ME-1633 on beamline ID19 at the European Synchrotron Radiation Facility.
Funding
This work is supported by ETH Zürich and the European Synchrotron Radiation Facility.
Declaration of interests
The authors report no conflict of interest.
Author contributions
M.C., G.S. and O.S. conceived the study and secured beam time for the X-ray experiments. M.C., L.P., G.S., B.L. and O.S. defined the experimental methodology. M.C., L.P., G.S., B.L., A.P. and O.S. performed the experiments. M.C. and L.P. post-processed and analysed the data. M.C. carried out the theoretical modelling. T.R. and D.W.M. performed the BEM simulations. M.C. wrote the initial draft of the paper. O.S., B.L. and A.R. provided resources. O.S. supervised the research. All authors contributed to the critical review and revision of the manuscript.
Data availability statement
The data that support the findings of this study are available upon reasonable request.
Appendix A. Derivation of the variational wetting relation
A perturbation of the bubble surface
$\delta {\boldsymbol{R}}$
induces a change in the wetting condition:

where
$\delta _{\perp } {\boldsymbol{n}}$
and
$\delta _{\parallel } {\boldsymbol{n}}$
represent the variations due to the normal
$\eta$
and tangential
$\tau$
components of the surface deformation, respectively. The term
$\delta \alpha$
denotes the change in the contact angle relative to the equilibrium state. Since the substrate is flat,
$\delta {\boldsymbol{n}}_{\textrm { s}} = 0$
. The normal unit vector is given by

where
$ \partial _{\theta } {\boldsymbol{R}} = {\partial {\boldsymbol{R}}}/{\partial \theta }$
and
$ \partial _{\phi } {\boldsymbol{R}} = {\partial {\boldsymbol{R}}}/{\partial \phi }$
. Considering the normal component of the surface deformation
$\delta _{\perp } {\boldsymbol{R}} = \eta {\boldsymbol{n}}$
, we obtain for the variation
$\delta _{\perp } {\boldsymbol{n}}$

The first term on the right-hand side is a vector orthogonal to
$\boldsymbol{n}$
, while the sum of the remaining two terms constitutes a vector parallel to
$\boldsymbol{n}$
. Given that the variation of a unit vector is orthogonal to the vector itself, the sum of the last two terms must be zero. Furthermore, since the angular coordinates
$\theta$
and
$\phi$
correspond numerically to the lengths of their respective arcs, it follows that
$\| \partial _{\theta } {\boldsymbol{R}} \times \partial _{\phi } {\boldsymbol{R}} \| = 1$
. Therefore,

Considering the tangential component of the surface deformation
$\delta _{\parallel } {\boldsymbol{R}} = \tau {\boldsymbol{t}}$
, we derive the variation
$\delta _{\parallel } {\boldsymbol{n}}$
by leveraging the directional derivative of the bubble surface:

Using the Frenet–Serret formulas,
$ \partial _{\theta } {\boldsymbol{n}} = - k{\boldsymbol{t}}$
, where
$k = -\|{\boldsymbol{R}}_0\|^{-1} = -1$
is the curvature of the bubble surface along the tangential direction
$\boldsymbol{t}$
, we obtain


Figure 20. First seven shape modes
$P_{l_{k}}^m(\cos (\theta ))$
for
$m=0$
describing the normal bubble deformation
$\eta$
for (a,b) a free bubble, (c,d) a pinned bubble with a static contact angle
$\alpha _0 = 3\pi /4$
and (e,f) a fixed-contact-angle bubble with
$\alpha _0 = 3\pi /4$
.
By inserting (A4) and (A6) in (A1) and considering that
${\boldsymbol{t}} \,{\boldsymbol\cdot}\, {\boldsymbol{n}}_{\textrm { s}} = - \sin (\alpha _0)$
,
${\boldsymbol{b}} \,{\boldsymbol\cdot}\, {\boldsymbol{n}}_{\textrm { s}} = 0$
and
$\tau |_{\theta =\alpha _0} = \cot (\alpha _0) \eta |_{\theta =\alpha _0}$
, we obtain the following variational wetting relation (3.4):

Appendix B. Visualisation of zonal shape modes for different bubble boundary conditions
Figure 20 presents the first seven shape modes
$ P_{l_k}^m(\cos (\theta ))$
for
$ m = 0$
, which describe the normal deformation
$ \eta$
of a bubble surface. The integer
$ k$
sequentially orders the permissible modes at a given
$ m$
. Figure 20(a,b) shows the shape modes for a free bubble. In this case, to comply with 2
$\pi$
-periodicity of
$\eta$
and its directional derivatives, both
$ l$
and
$ m$
must be integers. Consequently,
$l$
corresponds to
$k$
. Figure 20(c,d) illustrates the shape modes for a pinned bubble with a static contact angle
$ \alpha _0 = {3\pi }/{4}$
. Here, the surface deformation must vanish at the angle
$ \theta = \alpha _0$
. This boundary condition requires
$ l$
to be typically non-integer, meaning it no longer corresponds to
$ k$
. The zero-deformation line is represented by a bold black line. Figure 20(e,f) depicts the shape modes for a bubble with a fixed contact angle of
$ \alpha _0 = {3\pi }/{4}$
. In this scenario, the ratio between the surface deformation and its directional derivative along
$ \theta$
at
$ \theta = \alpha _0$
must equal
$ \cot (\alpha _0)$
. Like the pinned case, this condition typically results in non-integer values of
$ l$
, which again do not match
$ k$
. The
$ \cot (\alpha _0)$
line is shown as a bold black line.
Appendix C. Permissible values of shape-mode degrees for different bubble boundary conditions
Table 1 lists the permissible values of the degree
$ l$
of a shape mode for given indices
$ k$
and orders
$ m$
, in two different cases. Case (a) corresponds to a pinned wall-attached bubble with a static contact angle
$ \alpha _0 = {3\pi }/{4}$
, while case (b) corresponds to a fixed-contact-angle wall-attached bubble with
$ \alpha _0 = 3\pi /4$
. The values shown for case (a) correspond to the shape-mode spectrum depicted in figure 5(b), and those for case (b) correspond to the shape-mode spectrum depicted in figure 5(c).
Table 1. Permissible values of the degree
$l$
of a shape mode for a given set of indexes
$k$
and orders
$m$
in the cases of (a) a pinned wall-attached bubble with a static contact angle
$\alpha _0 = 3\pi /4$
and (b) a fixed-contact-angle wall-attached bubble with
$\alpha _0 = 3\pi /4$
.

Appendix D. Derivation of the pressure term accounting for the rigid boundary using the method of images
The radial motion of a spherical bubble of radius
$R(t)$
in an unbounded medium under spatially uniform acoustic forcing can be described using the Rayleigh–Plesset equation for mildly compressible Newtonian media, which reads

To incorporate the influence of the solid boundary, we employ the method of images (Strasberg Reference Strasberg1953). This approach, rooted in the potential flow theory, simulates a flat boundary within the flow field by introducing a virtual bubble mirroring the real bubble on the opposite side of the boundary. Due to the linear nature of the potential flow approximation, the influence of each bubble on the flow field can be analysed independently. Following the approach adopted in several other studies (Zabolotskaya Reference Zabolotskaya1984; Oguz & Prosperetti Reference Oguz and Prosperetti1990; Mettin et al. Reference Mettin, Akhatov, Parlitz, Ohl and Lauterborn1997; Pelekasis et al. Reference Pelekasis, Gaki, Doinikov and Tsamopoulos2004; Bremond et al. Reference Bremond, Arora, Ohl and Lohse2006), the velocity field
$u_{ { l}}(r,t)$
induced by each bubble can be derived from the continuity equation and is expressed as

where
$r$
is the radial coordinate, with the origin of the coordinate system at the centre of the bubble. The corresponding pressure field
$p_{ { l}}(r,t)$
can be derived from the momentum equation:

where the nonlinear term is omitted due to its relatively smaller magnitude compared with the other terms. Integrating this equation yields the pressure field:

which acts as an additional driving pressure for the second bubble. Therefore, the equation governing the radial motion of a bubble in proximity to a solid boundary can be approximated as follows (5.3):

Appendix E. Zhou model for the bubble gas pressure
The method is based on the assumption that the gas pressure within the bubble is uniform (Prosperetti, Crum & Commander Reference Prosperetti, Crum and Commander1988). For an ideal gas, the radial velocity
$u_{ { g}}(r)$
is expressed as

leading to an exact relation for the gas pressure:

where
$r$
is the radial coordinate,
$\gamma$
is the gas specific heat ratio,
$K_{{ g}}$
is the gas thermal conductivity and
$T_{{ g}}$
is the gas temperature. The temperature profile
$T_{{ g}}(r)$
inside the bubble is divided into three regions: (i) an internal layer with uniform temperature, (ii) a buffer layer and (iii) an external layer with a linear temperature distribution. The change in bubble surface temperature is negligible, and thus
$T_{{ g}}|_{R} \approx T_{ { l}}$
(Prosperetti et al. Reference Prosperetti, Crum and Commander1988). The volume-averaged temperature
$T_{ { g}_{\textit {i}}}$
in each region
$i$
is computed using the ideal gas law:

where
$\rho _{{g}_{i}}$
is the volume-averaged gas density in region
$i$
and
${\mathcal{R}}$
is the specific gas constant. The gas density
$\rho _{{ g}_{i}}$
can be computed using the continuity equation for each region as

where
$m_{ { g}_{\textit {i}}}$
is the gas mass in region
$i$
and
$f_{j}$
is the mass flux across the interface
$j$
:

with
$\rho _{ { g,}{j}}^{\textit{u}w}$
being the density on the upwind side of the interface
$j$
,
$u_{ { g,}{j}}^{\textit{rel}}$
the convective velocity (difference between real gas velocity and cell interface velocity) and
$S_j$
the interface surface area.
Appendix F. Implementation of the Francescutto & Nabergoj model for wall- attached bubbles
The instability pressure threshold
$p_{ {th}}, {l}$
for a shape mode with integer degree
$l$
, and consequently for a shape mode with order
$m=l$
, can be obtained under the assumption of linear oscillations by determining the stability condition of the Mathieu equation for the shape-mode amplitude and reads

where
$C_{{th},l}$
is the relative bubble oscillation amplitude:

in which the coefficients
$e$
,
$f$
and
$g$
are given by



The dissipative effects are taken into account by the term
$\delta$
, which reads

where
${\,\rm \mu} _{ { e}}$
represents the effective viscosity of the liquid, incorporating the combined contributions from viscous, thermal and acoustic radiation losses. Within the bubble radius range examined in this study (
${60}{-}{140}{\;\unicode{x03BC}} {\textrm{m}}$
) and under linear oscillation conditions, the effective viscosity is approximately 30 to 70 times greater than the contribution from liquid viscosity
${\,\rm \mu} _{ { l}}$
alone, with these values corresponding to the smallest and largest bubbles, respectively; all values, including intermediate ones, were extracted from the theoretical plots in Chapman & Plesset (Reference Chapman and Plesset1971). Additional viscous losses may arise due to the bubble proximity to the wall and the motion of the contact line; however, these effects are difficult to quantify. Despite neglecting them, the model shows good agreement with experimental data, suggesting that their contribution is relatively minor. Finally,
$\omega _{\textit{res}} = 2\pi f_{{ res}}$
is the resonance angular frequency of the bubble, and
$\omega _{ { d}} = 2 \pi f_{ { d}}$
is the ultrasound driving angular frequency. The bubble resonant frequency for linear oscillations is numerically determined by solving the fully nonlinear radial dynamics model (5.3), incorporating the Zhou gas model, to properly account for the thermal behaviour of the bubble. A small driving pressure (
$p_{ { a}}={1}\,{\textrm{kPa}}$
) is used to ensure the bubble remains in the linear oscillation regime. The resonant frequency is identified using a golden ratio search algorithm, which finds the driving frequency that maximises radial expansion for each bubble radius. The resonant bubble radius corresponding to the ultrasound driving frequency employed,
$f_{ { d}} = {30}\,\textrm{kHz}$
, is
$R_0 = {85.8}\;\unicode{x03BC} \textrm{m}$
.
Appendix G. Implementation of the BEM for wall-attached bubbles
Given the Mach number
$M=\dot {R}/c_{ { l}} \sim O(10^{-3})$
, Reynolds number
$Re=2\rho _{ { l}}\dot {R}R/{\,\rm \mu} _{ { l}} \sim O(10^{2})$
and Bond number
$Bo = (\rho _{ { l}} - \rho _{ { g}}) g R^2/\sigma \sim O(10^{-3})$
for the problem at hand, we can describe the flow velocity
$\boldsymbol{u}_{ { l}}(\boldsymbol{x},t)$
and pressure
$p_{ { l}}(\boldsymbol{x},t)$
within the liquid domain, denoted by
$\Omega$
, using the incompressible Euler equations without considering gravitational forces:

Under these conditions, the flow outside the boundary layer can be approximated as irrotational, i.e.
$\boldsymbol{\nabla } \times \boldsymbol{u}_{ { l}}= 0$
. Since the domain
$\Omega$
is simply connected, we can introduce a scalar potential
$\phi (\boldsymbol{x},t)$
such that
$\boldsymbol{\nabla }\phi = \boldsymbol{u}_{ { l}}$
. This simplifies Euler equations into the Laplace equation for the potential and the Bernoulli equation:

By leveraging the Green’s third identity (derived from Green’s theorem), which relates the values of a harmonic function (such as the potential
$\phi$
) at any point inside the domain to the values of the function and its normal derivative on the boundary, we can express the potential
$\phi$
at any point
$\boldsymbol{x}$
within the domain in terms of an integral over the liquid–gas interface
$S$
, as follows:

where
$c(\boldsymbol{x})$
is the solid angle and the Green’s function is given by
$g(\boldsymbol{x},\boldsymbol{y}) = 1/||\boldsymbol{x} - \boldsymbol{y}||$
(e.g. Steinbach Reference Steinbach2008). Thus, the dimensionality of the problem is reduced from solving for values of the potential in the entire domain to solving only for the values along the boundary. To solve numerically the boundary integral relation, we discretise the bubble surface using a triangular mesh, where the potential values
$\phi (\boldsymbol{x}_i)$
and the normal gradients
$\partial \phi (\boldsymbol{x}_i)/\partial {\boldsymbol{n}}$
at the triangle vertices
$\boldsymbol{x}_i$
serve as the numerical degrees of freedom. For computing the surface integral
$\int _S\ldots \mathrm{d}S(\boldsymbol{y})$
over this triangular grid, we utilise a seven-point two-dimensional quadrature rule for each triangle (Radon Reference Radon1948, p. 295). To obtain the values of
$\phi (\boldsymbol{y})$
and its normal derivative
$\partial \phi (\boldsymbol{y})/\partial {\boldsymbol{n}}$
at surface points
$\boldsymbol{y}$
within a triangle, bilinear interpolation is employed. To address the singularities that occur in the evaluation of
$g(\boldsymbol{x}_i,\boldsymbol{y})$
when the surface point
$\boldsymbol{y}$
coincides with the collocation point
$\boldsymbol{x}_i$
, we use a polar coordinate transform on a unit triangle (Kieser, Schwab & Wendland Reference Kieser, Schwab and Wendland1992) in combination with a seven-point one-dimensional Gauss–Legendre quadrature rule. We calculate the solid angle
$c(\boldsymbol{x})$
at the collocation point
$\boldsymbol{x}_i$
by applying the
$4\pi$
rule, as described in Klaseboer et al. (Reference Klaseboer, Fernandez and Khoo2009, equation (6)). Finally, we model the rigid wall as a mirror bubble, in accordance with the potential flow theory.
At the contact line
$\boldsymbol{\Gamma }$
between the bubble surface and the solid wall, a boundary condition is imposed that allows the contact line to move along the wall plane with a velocity
$w$
, which is determined by the difference between the cosines of the equilibrium contact angle
$\alpha _0$
and the instantaneous contact angle
$\alpha (t)$
. This relationship, derived from Ren & E (Reference Ren and Weinan2007, equation (28)), is described by

where
$\beta$
is the effective friction coefficient. This formulation leads to a condition for the normal derivative of the scalar potential
$\partial \phi (\boldsymbol{x}) / \partial {\boldsymbol{n}}$
along the contact line. To ensure numerical stability, we prevent the angle between the bubble surface and the wall from becoming too small – which corresponds to
$\cos (\alpha (t))$
approaching either
$1$
or
$-1$
– by introducing an additional nonlinear term:

This term is designed as the simplest rational function of
$\cos (\theta )$
that diverges to
$+\infty$
as
$\cos (\theta ) \to -1$
and to
$-\infty$
as
$\cos (\theta ) \to +1$
, is antisymmetric about the origin and has a zero slope at
$\cos (\theta ) = 0$
. The
${\sigma _{{ lg}}}/{\beta }$
ratio is calibrated to 0.2 to replicate the contact-line mobility observed in experiments. The parameter
$c_n$
is set to 0.01, the minimum value ensuring stable simulations. These parameters are maintained constant in all simulations.
Ultimately, we obtain a dense linear system of equations whose unknowns are the normal gradients
$\partial \phi (\boldsymbol{x})/\partial \boldsymbol{n}|_{\boldsymbol{x}=\boldsymbol{x}_i}$
at vertices away from the contact line and the potential values
$\phi (\boldsymbol{x}_i)$
at vertices on the contact line. We solve this linear system using the BiCGStab solver (van der Vorst Reference van der Vorst1992). By combining the normal potential gradients at the interface
$\partial \phi (\boldsymbol{x})/\partial \boldsymbol{n}|_{\boldsymbol{x}=\boldsymbol{x}_i}$
with the tangential gradients at the interface, which can be computed from
$\phi (\boldsymbol{x}_i)$
, we can numerically estimate the advection velocity at the interface nodes,
$\boldsymbol{\nabla } \phi (\boldsymbol{x}_{i}) = u_{\mathrm{l}}(\boldsymbol{x}_{i})$
. Together with the contact line velocity
$w$
, this enables us to evolve the interface over time by advecting the nodes. We evolve the potential on the bubble surface over time using the Bernoulli equation (e.g. Blake, Tomita & Tong Reference Blake, Tomita and Tong1998, equation (8)):

where
${\textrm D} /{\textrm D}t = \partial /\partial t + \boldsymbol{u}_{ { l}}\,{\boldsymbol\cdot}\, \boldsymbol{\nabla }$
is the material derivative and
$\kappa$
is the local bubble curvature, which is computed by using the scheme of Rusinkiewicz (Reference Rusinkiewicz2004, § 3.1). To reduce computational complexity, we assume that the gas inside the bubble undergoes an adiabatic process. Although this model may not accurately represent the breathing mode, as demonstrated in § 5.1, its influence on the evolution of shape modes appears to be minimal, as shown in the following. Here,
$V(t)$
represents the bubble volume at time
$t$
, and
$V_0$
denotes the initial bubble volume.
To accurately evolve grid points and their associated potentials over time, we apply a four-stage Runge–Kutta scheme coupled with an adaptive time-stepping discretisation method. Each simulation begins with a cut geodesic icosahedron of radius
$R_0$
, positioned at a height of
$-R_0\cos {\alpha }$
below the wall, with all initial potential values set to zero. The triangle mesh typically comprises several thousand elements. After a number of time steps, the mesh quality is recalibrated using standard methods described in § 6.5.3 of Botsch et al. (Reference Botsch, Kobbelt, Pauly, Alliez and Levy2010). These methods include vertex relaxation, edge flipping and the splitting or collapsing of edges based on a specified target edge length. In our approach, this target length is determined locally, taking into account both the local curvature and potential gradient. To prevent abrupt changes in the target edge length between time steps, we apply a moving time average for the target length.