Ultrasound-induced nonlinear oscillations of a spherical bubble in a gelatin gel

Abstract We experimentally observe and theoretically analyse oscillations of a spherical bubble in a gelatin gel under ultrasound irradiation to quantify viscoelastic effects on the nonlinear bubble dynamics. A bubble nucleus is generated by focusing a laser pulse into a 6 wt% gelatin gel supersaturated with dissolved air, which enables us to control the bubble radius at mechanical equilibrium via influx of the gas air into the bubble. Linearized and finite-amplitude oscillations of the bubble are driven by 28 kHz ultrasound and recorded by a high-speed camera; the resonance curves of the oscillation amplitude as a function of the equilibrium radius are constructed for different ultrasound intensities. First, the viscosity and shear modulus of the gel are obtained by fitting the resonance curve (for the lowest ultrasound intensity) to the linearized solution of the Rayleigh–Plesset model that accounts for the gel's nonlinear elasticity of neo-Hookean type and diffusive effects on the bubble dynamics. Next, finite-amplitude oscillations of the bubble are compared with the nonlinear Rayleigh–Plesset calculations. The comparison suggests a need to include the gel's elasticity in the calculations to more accurately reproduce the nonlinear bubble dynamics. Another important finding is that the so-called spring softening feature appears in the experimentally determined resonance curve as the oscillation amplitude increases, which can be predicted by the Rayleigh–Plesset model. Furthermore, our experiment with the highest ultrasound intensity shows non-spherical oscillation of mode 1 that does not appear in the case of water but can be predicted by shape instability theory.


Introduction
Acoustically driven bubbles play an important role in therapeutic ultrasound applications such as histotripsy (Coussios & Roy 2008;Brennen 2015;Bader, Vlaisavljevich & Maxwell 2019). Large-amplitude volume oscillations of bubbles under high-intensity megasonic wave irradiation often result in extremely high pressure and temperature of the bubble contents at violent collapse, which can destroy malignant tumours. Since human tissues are viscoelastic, a better understanding of the role of viscoelasticity in acoustic bubble dynamics is desirable. Although a number of theoretical and numerical studies on bubble dynamics in soft matter have been reported (Brujan 2011;Dollet, Marmottant & Garbin 2019), experimental data are yet lacking, and thus limiting model validation as well as quantitative evaluation of the viscoelasticity effects on bubble dynamics. In particular, acoustically induced nonlinear oscillations of a bubble in soft matter have not been systematically observed, and the corresponding nonlinear bubble dynamics remain elusive. For this reason, we experimentally study the nonlinear oscillations of a bubble in a gelatin gel under ultrasound irradiation and examine the role of viscoelasticity.
Essential features of bubble dynamics can be extracted from Rayleigh-Plesset-type equations under the assumption that bubbles keep their spherical shape (Leighton 1994;Brennen 1995). The effect of viscoelasticity of the medium surrounding bubbles can be incorporated into the Rayleigh-Plesset modelling (Prosperetti 1982;Yang & Church 2005;Gaudron, Warnez & Johnsen 2015;. Yang & Church (2005) studied ultrasound-induced nonlinear dynamics of spherical bubbles in linear Kelvin-Voigt solids and showed that the elasticity can suppress the oscillation amplitude.  examined more carefully effects of the viscoelasticity including stress relaxation on bubble dynamics. Gaudron et al. (2015) coupled the Keller-Miksis equation, describing dynamics of spherical bubbles in compressible media (Keller & Miksis 1980;Prosperetti & Lezzi 1986), to a Kelvin-Voigt model with (nonlinear) neo-Hookean elasticity and showed that nonlinear elasticity can lead to milder bubble collapse.
Experimental observations of spherical cavitation bubbles in soft matter are required to validate the Rayleigh-Plesset-type equations equipped with constitutive equations of viscoelastic solids. Since acoustic bubble nucleation is a random event (Oguri & Ando 2018), the main target in most of the previous experiments is bubbles that are nucleated by focusing a laser pulse into a specific location (see, for example, Brujan & Vogel (2002)). With laser-nucleated bubbles in soft matter, one can observe both free and forced bubble oscillations with the aid of a high-speed camera. Estrada et al. (2018) observed free oscillations of laser-induced bubbles in a polyacrylamide gel and compared them with Rayleigh-Plesset calculations with a Kelvin-Voigt viscoelastic model; the nonlinear oscillations up to the third rebound are found to be reproduced by the model. Hamaguchi & Ando (2015) observed small-amplitude oscillations of spherical bubbles in a gelatin gel under 28 kHz ultrasound irradiation and compared them with the linearized solution of the Rayleigh-Plesset equation with the linear Voigt model; this comparison allows one to obtain viscosity and elasticity coefficients of the gel under the ultrasound forcing. Experimental observations of finite-amplitude nonlinear bubble oscillations under ultrasound irradiation are essential in the context of the therapeutic ultrasound, but yet are lacking.
The goal of this study is to examine the role of viscoelasticity in finite-amplitude oscillations of spherical bubbles in a viscoelastic solid under ultrasound irradiation. While a megasonic wave is often used in therapeutic ultrasound, we use a lower ultrasound frequency (28 kHz), which enables us to obtain sufficient temporal and spatial resolutions of the bubble dynamics. Following the approach of Hamaguchi & Ando (2015), we experimentally observe both small-and finite-amplitude oscillations of laser-nucleated bubbles in a gelatin gel under ultrasound irradiation. We use the small-amplitude bubble oscillations to obtain the gel's viscoelastic properties from their comparison with the linearized Rayleigh-Plesset solution; ultrasound intensity and the radius of bubbles at mechanical equilibrium are treated as parameters, respectively, to examine the extent of nonlinearity and resonance. We also examine the finite-amplitude bubble oscillations (with varying ultrasound intensity and equilibrium bubble radii) and compare them with Rayleigh-Plesset calculations that account for nonlinear elasticity of the gel (Gaudron et al. 2015). To the authors' knowledge, this is the first time a nonlinear dynamics model for spherical bubbles in gels has been validated by a comparison with experiments.
In what follows, we first introduce theoretical and experimental methods to study spherical bubble dynamics in viscoelastic solids and then discuss the role of viscoelasticity through comparison between the theory and experiments.

Full model
A model for spherical, nonlinear bubble dynamics in viscoelastic solids can be derived based on continuum mechanics (Keller & Miksis 1980;Prosperetti, Crum & Commander 1988;Gaudron et al. 2015;Estrada et al. 2018;Mancia et al. 2021). Oscillations of a spherical bubble in terms of instantaneous radius R(t) can be described by the Keller-Miksis equation where liquid compressibility in the far-field is taken into account (Keller & Miksis 1980;Estrada et al. 2018), where ρ 0 and c 0 are the (undisturbed) density and sonic speed of the medium surrounding the bubble, respectively, the dot represents a temporal derivative, and σ is the surface tension. The far-field pressure consists of atmospheric pressure and acoustic pressure p ∞ (t) = p 0 + p a (t). The stress integral S for the finite-deformation Kelvin-Voigt model is obtained as (Gaudron et al. 2015;Hamaguchi & Ando 2015;Estrada et al. 2018) where μ and G are the shear viscosity and shear modulus of the surrounding medium, respectively, and l is the bubble radius at strain-free state. We assume that the bubble consists of non-condensable gases (say, air) only (with a negligible amount of vapour). The homobaric bubble pressure p g (t) obeys (Prosperetti et al. 1988 where γ is the ratio of specific heats (γ = 1.4 for air), K(T) is the thermal conductivity of the gas, and the radial coordinate r is measured from the bubble centre. We use the temperature-dependent thermal conductivity of air (Prosperetti et al. 1988), The temperature field T(r, t) inside the bubble is determined by solving where C p is the specific heat of the gas at constant pressure. The velocity field inside the bubble u g (r, t) is given by (2.6) Equations (2.3) to (2.6) are valid when the gas inside the bubble follows the ideal gas law p g = ρ g R g T, where R g is the gas constant. We employ the cold liquid assumption (Prosperetti et al. 1988;Preston, Colonius & Brennen 2007), so that the liquid temperature is undisturbed at room temperature T ∞ (i.e. T(R, t) = T ∞ ). In addition, the insulated boundary condition is applied at the bubble centre (i.e. ∂T/∂r(0, t) = 0). The bubble pressure at mechanical equilibrium (but with prestrain in the surrounding medium) is where R 0 is the bubble radius at mechanical equilibrium. Note that the equilibrium radius l at which the medium surrounding the bubble is strain free allows one to account for prestrain effects on the bubble pressure. In the present experiment, the bubble nuclei are generated by focusing a pulsed laser into (originally strain free) gels of 6 wt% gelatin concentration (see § 3.1) so that prestrain will be built up after the laser focusing; ratio on the strain-free radius to the equilibrium radius is set to l/R 0 = 0.6, which is close to the prediction of Ando & Shirota (2019). In summary, the nonlinear oscillations of a spherical bubble in a viscoelastic medium are described by solving (2.1), (2.3) and (2.5) simultaneously. The model is solved numerically: we semidiscretize the partial differential equation (2.5) by central differences (Prosperetti et al. 1988), and integrate in time using the Cash-Karp Runge-Kutta method (Press et al. 2007). We use a computational domain y = r/R(t) where the number of a uniform spatial mesh is set to 200 with which the numerical solutions converge enough in our problems. The initial conditions are given by R(0) = R 0 ,Ṙ(0) = 0, p g (0) = p g,0 and T(r, 0) = T ∞ . The physical properties used in this model are summarized in table 1; the value of surface tension σ is taken from Kikuchi, Sakai & Takagi (1991), and the viscosity μ and shear modulus G are fitted values (see § 4.1).

Linearized solution
The linearized solution to the full model equations (2.1) to (2.7) can be derived by perturbation techniques, which enables one to obtain the viscosity and shear modulus of the surrounding medium from its comparison with small-amplitude bubble oscillation data (Hamaguchi & Ando 2015). Given sinusoidal acoustic pressure p a (t) = C cos(2πft) with sufficiently small amplitude C and frequency f , the bubble response is expected to be linearized (R = R 0 (1 + x) where |x| 1) and the bubble content tends to  behave polytropically. In this limit, the full model equations reduce to the linear ordinary equation of a mass-damper-spring system where mass, damping and spring constants are given, respectively, by Here, κ eff and μ eff are the effective polytropic index and viscosity (Chapman & Plesset 1971;Prosperetti 1977a;Ando, Colonius & Brennen 2009), where Υ is the complex transfer function of the Péclet number Pe = R 2 0 (2πf )/α g , Note that the nonlinear elasticity term of prestrain still remains in (2.7) and (2.9a-c), regardless of whether the bubble response is linear or nonlinear. When the prestrain is assumed to vanish (i.e. l = R 0 ), the elastic model reduces to the linear model of Hamaguchi & Ando (2015). The steady-state solution of the linearized equation (2.8) is given by where X and δ are the oscillation amplitude and phase shift, respectively, (2.14) Here, f N and ζ are the natural frequency and damping ratio, respectively, with which the resonant frequency f R can be defined as In the present experiment, the equilibrium radius of bubbles, R 0 , is varied but the ultrasound frequency f is fixed at 28 kHz. When it comes to constructing resonance curves of bubble oscillations, we can plot the oscillation amplitude X as a function of the equilibrium bubble radius R 0 , for the natural frequency f N and equilibrium radius R 0 are one-to-one through (2.9a-c) and (2.15a,b).

Generation of a spherical gas bubble nucleus in a gelatin gel
We use a 6 wt% air-supersaturated gelatin gel as a viscoelastic solid. A sol generated by mixing powdered gelatin (G 2500, Sigma-Aldrich) and heated deionized water is poured into an acrylic case (38 mm × 38 mm × 15 mm), whose bottom is sealed with a thin plastic wrap. The sol turns to a gel after overnight cooling at 4 • C. After returning to room temperature, the gel is supersaturated with dissolved air (Hamaguchi & Ando 2015) since mass diffusion is much slower than heat diffusion. While viscoelastic properties of gels under strain rates up to approximately 100 Hz are measured commonly by commercial rheometers, those of the gel deformed under ultrasound irradiation cannot be measured in the same manner (Kwon & Subhash 2010;Jamburidze et al. 2017;Saint-Michel & Garbin 2020) and thus need to be treated as unknown parameters in this study.
A spherical bubble nucleus can be generated by focusing a laser pulse into the gelatin gel as illustrated in figure 1(a). We use a Q-switched Nd:YAG laser (ULTRA 50 GRM, Quantel) delivering an infrared laser pulse whose wavelength and duration are 1064 nm and 6 ns, respectively. The laser pulse is focused into the gel through a microscope objective (×40, NA = 0.6). A micron-sized cavity with dissolved air appears as a result of laser-induced plasma and subsequent bubble formation. Since the surrounding gel is supersaturated with the dissolved gas, the cavity exhibits gradual growth as a spherical gas bubble as a result of the gas influx, see figure 1(b). The bubble growth is obtained approximately over one hour as shown in figure 1(c) and is thus so slow that (unsteady) dynamics of the bubble are totally unimportant. Namely, we can control the radius of the bubble at mechanical equilibrium, R 0 , based on the mass diffusion.
It is instructive to report on the selection of gelatin concentration at 6 wt%. To confirm sphericity of the laser-induced gas bubble nuclei, we recorded the side view of the bubble in the gels of the 6 wt% and 12 wt% gelatin concentration, see figure 2. In the case of the higher gelatin concentration, deviation of the bubble shape from a sphere is further emphasized, given the same objective lens with NA = 0.6. Here, we define the bubble's roundness as 4π × area/( perimeter) 2 where 'area' and 'perimeter' are calculated from the (two-dimensional) images of the bubbles (Hamaguchi & Ando 2015). The roundness of the bubbles in the 6 wt% gel is 1.05 ± 0.01 (out of five measurements including figure 2a), while that in the 12 wt% is 1.11 ± 0.06 (out of five measurements including figure 2(b). For the bubble in the stiffer gel to be spherical, there is a need to use objective lenses with higher NA that enable pointwise plasma formation. However, the  working distance of higher-NA objective lenses is shorter, making it more difficult to create bubbles away from the container's wall. Since the observation target of our concern is isolated bubble dynamics, laser-induced bubble nuclei need to be created sufficiently away from any boundaries. This technical difficulty hinders the study of isolated bubble dynamics in stiffer gels with higher gelatin concentration. The case with the 3 wt% gelatin concentration, on the contrary, will lead to more spherical bubble nuclei but give rise to shape instability more easily at smaller-amplitude oscillations (Hamaguchi & Ando 2015). Thus, it is also an experimental challenge to obtain spherical oscillations of finite amplitude in the lower concentration case.

Recording bubble oscillations driven by ultrasound
The bubble oscillations driven by underwater ultrasound from a 28 kHz planer piezoelectric transducer (0.028Z45I, JAPAN PROBE, wavelength in water λ = 54 mm) are recorded by a high-speed camera (FASTCAM SA-X2, Photron, frame rate 300 kfps, resolution 2.88 μm pixel −1 , exposure time 293 ns) with backlighting of LED illumination 'Air-supersaturated' gel with a spherical bubble LED illumination High-speed camera (300 kfps, 2.88 µm pixel -1 , exposure time: 293 ns) Piezoelectric transducer (28 kHz, wavelength in water λ = 54 mm) (SLG-150V, REVOX) and a zoom lens (×300, CX-5040SZ, HIROX), as illustrated in figure 3. The transducer is attached on the bottom of an acrylic container filled with tap water. The distance between the transducer and the free surface of the water is set to 5λ/4 in order to form a standing pressure wave (Yamashita & Ando 2020). The bubble generated in the gel (see § 3.1) is located at the antinode of the standing wave (i.e. λ/2 above the transducer). Bitmap images taken from the high-speed recordings are converted to binary data using the graythresh function in MATLAB and the area of the projected bubble is obtained by the regionprops function to determine the area-equivalent radius.
There are two important parameters in our experiments: (i) radius of the bubble at mechanical equilibrium R 0 (see § 3.1); (ii) ultrasound pressure amplitude p A .
The acoustic pressure p a (t) in the water (without the gel sample) is measured by a hydrophone (HCT0310, Onda). The acoustic pressure in the gel is expected to be similar to the underwater pressure, for their acoustic impedances are similar. In the experiment, the pressure amplitude p A = max( p a ) is set either at 0.03 atm (for linearized oscillation) or at 0.24 atm (for finite-amplitude oscillation).
In the high-speed camera recording, we consistently use the same gas bubble nucleus created in a specific position of the gel sample. Note that the structure of gels is inhomogeneous in space and the gel viscoelasticity field is essentially non-uniform (Shibayama 1998). This approach thus enables us to eliminate effects of the gel's inhomogeneity on the analysis of bubble dynamics for each value of R 0 .

Small-amplitude bubble oscillations
First, we consider the case of small-amplitude bubble oscillations under low-amplitude ultrasound pressure p A = 0.03 atm. In figure 4, as a representative example, we present the result for the bubble (under resonance f /f R = 1.00) whose equilibrium radius is dimensionally at R 0 = 115.9 μm. It follows from a series of recorded images (figure 4a) that the bubble keeps its spherical shape for the duration of our observation. From the evolution of the area-equivalent bubble radius (figure 4b), it turns out that the oscillation amplitude in the radius during both the growth and shrinkage phases is approximately 0.05R 0 . We may thus say that the oscillation amplitude is sufficiently small for the dynamics to be linearized.
We repeated the same experiment as in figure 4, but varying the equilibrium radius R 0 of the same gas bubble nucleus based on the diffusion-driven growth ( § 3.1), while the ultrasound frequency and amplitude were both fixed. This procedure allows for constructing the resonance curve, i.e. the relation between the oscillation amplitude max(R) and the equilibrium bubble radius R 0 as in figure 5 where the distinct experimental results with different R 0 are denoted by circles. Now that the oscillation amplitude is small and the dynamics are expected to be linearized as supported by figure 4 (for the case of the resonant bubble), we compare the resonance curve with linear theory ( § 2.2). In this comparison, we use MATLAB's built-in function lsqcurvefit for unknown shear viscosity μ and shear modulus G of the gel to be fitted. The fitted solution is also shown in figure 5, which reasonably agrees with the experimentally determined resonance curve. The fitted values of the shear viscosity and shear modulus are μ = 18.3 mPa s and G = 4.0 kPa, respectively, and used for analysis in the succeeding sections. Finally, we note that these fitted values are lower than those reported in Hamaguchi & Ando (2015), for viscosity and rigidity coefficients of the gel are sensitive to its production and storage processes (te Nijenhuis 1981) and can vary between different gel samples even at the same gelatin concentrations.

Finite-amplitude bubble oscillations
Next, we consider the case of finite-amplitude bubble oscillations under the high ultrasound pressure p A = 0.24 atm. In figure 6(a), we present a visualization example of spherical oscillation of the bubble with R 0 = 124.0 μm (and f /f R = 1.06). The evolution of the corresponding area-equivalent bubble radius is shown in figure 6(b), illustrating that the amplitude during the growth phase is approximately 0.3R 0 but that during the shrinkage phase is smaller, which is a sign of the nonlinearity. The full model predictions ( § 2.1) are also shown in this figure. In the simulation, we use the shear viscosity and shear modulus obtained from the linear analysis ( § 4.1); for comparison, we also consider the case of no elasticity (G = 0). It follows from the comparison that the simulation agrees well with the experiment but that with the assumption of no elasticity gives rise to an underestimate in the oscillation amplitude; we discuss this trend in § 4.3. In short, the nonlinear dynamics of the spherical bubble in this particular example can be reproduced by the Rayleigh-Plesset-type model that is coupled to a Kelvin-Voigt constitutive equation with neo-Hookean elasticity, provided values of the shear viscosity and shear modulus are properly tuned, say, based on linear analysis. It is instructive to report on the larger-amplitude bubble oscillations even under the same pressure amplitude p A = 0.24 atm. In figure 7(a), we present a visualization example of the bubble oscillations with R 0 = 107.8 μm (and f /f R = 0.93). In this case, the bubble initially shows spherical oscillations (with larger amplitude in comparison with figure 6), but later non-spherical oscillations of mode 1 appear (in addition to the volume oscillations of mode 0) after t * ≈ 13.0. Note that non-spherical oscillations of mode 1 cannot happen in water with no rigidity (Prosperetti 1977b;Versluis et al. 2010;Liu et al. 2012), but can appear in our case, for the gelatin gel possesses rigidity that works as restoring force in the translational direction, as predicted by the analysis of Murakami, Gaudron & Johnsen (2020). In figure 7(b), we present the evolution of the corresponding area-equivalent bubble radius and compare it with the full model predictions. As in figure 6(b), we obtain reasonable agreement between the experiment and simulation, provided the bubble oscillates spherically. Once non-spherical oscillations appear after t * ≈ 13.0, the oscillation amplitude in volume is reduced, for some of the inputted ultrasound energy is used in the shape oscillations (Feng & Leal 1997). In this case, the Rayleigh-Plesset simulation for spherical bubble dynamics gives rise to an overestimate in the volume oscillation amplitude.
To understand the origin of the nonlinear response, we performed supplemental simulation for the same conditions as in figure 7 but with the linear elastic model (Hamaguchi & Ando 2015); see figure 9 in Appendix A. The simulation result with in larger-amplitude oscillations of spherical bubbles, particularly at radii smaller than equilibrium (Gaudron et al. 2015). For these effects to be visible, a far larger far-field pressure forcing (say, at tenfold amplitude p A in comparison with figure 7) would be necessary; see figure 10 in Appendix A. However, such large pressure forcing would give rise to significant departures from spherical oscillations. We repeated the same experiment with the pressure amplitude fixed at p A = 0.24 atm but with varying the equilibrium bubble radius R 0 , which allows for constructing the resonance curve (including the data from figures 6 and 7) for the finite-amplitude bubble oscillations as presented in figure 8. For easier comparison, the resonance curve of the small-amplitude oscillations is also included in this figure. The resonance curve of the large-amplitude oscillations shows nonlinear characteristics. To be specific, the peak is shifted to the smaller equilibrium radius or the resonant bubble radius is reduced in other words. This nonlinear feature is known as spring softening and is reported in numerical studies for bubbles in water (Lauterborn 1976;Fujiwara & Shima 1980;Kubo et al. 2012). While bubbles oscillating in water give rise to fission (Brennen 2002), the bubbles in the gel are rather stable against fission, allowing for experimental confirmation of the spring softening feature.

The role of viscoelasticity
Finally, we conduct parametric studies on the resonance curve with the shear viscosity and shear modulus to discuss the role of viscoelasticity. In figure 8, the solid blue line is the resonance curve obtained from the full model with the gel's shear viscosity and shear modulus μ = 18.3 mPa s and G = 4.0 kPa. In addition, the resonance curves obtained from the full model with no shear viscosity μ = 0 and doubled shear viscosity μ = 36.6 mPa s (but with the same shear modulus) are shown as red lines, and those with no shear modulus G = 0 and doubled shear modulus G = 8.0 kPa (but with the same shear viscosity) are shown as green lines. Shear viscosity suppresses the oscillation amplitude, thus reducing the soft spring feature and increasing the resonant bubble radius. Shear modulus increases the resonant bubble radius, while it has a small impact on the oscillation amplitude.
The discrepancy in the resonance curves for the higher-intensity ultrasound between the experiment and the full model (μ = 18.3 mPa and G = 4.0 kPa) is found to be emphasized especially for the case of the bubble whose equilibrium radius is smaller than the resonant radius. We note that all the experimental data reported in figure 8 are based on the same bubble nucleus and the order of the experimental runs is from the smallest to largest size of the bubble growing under gas supersaturation, with alternately using the higher-and lower-intensity ultrasound. We thus speculate that such a trend in the discrepancy is due to reduced viscosity and rigidity of the gel structure damaged by the larger bubble oscillations at the resonant condition (Movahed et al. 2016). For example, the gel in figure 4 (corresponding to the red circle in figure 8) had already been damaged by larger-amplitude bubble oscillations at the stage of figure 7 (corresponding to the red filled circle in figure 8). This implies that our fitted viscoelastic parameters dominantly correspond to the damaged gel. This supports the observation that the full model with the higher viscoelastic coefficients (chain lines) yields better agreement with the experimental data for the smaller bubble near which the gel was not yet damaged significantly.
At this stage, we revisit the simulation results in figures 6(b) and 7(b), respectively, for the case of smaller-and larger-amplitude oscillations. In these simulations, we consistently use a viscosity of μ = 18.3 mPa s (fitted from the resonance curve in figure 5). However, there is no need to stick to the fitted value, for the gel damage is expected to occur during the repeated experiments with the same bubble oscillating at finite amplitude. To quantitatively illustrate the viscous effect on the bubble dynamics, we perform supplemental simulations for the cases of μ = 0 (inviscid) and μ = 36.6 mPa s (doubled) and compare them with that with μ = 18.3 mPa s in figures 11 and 12, respectively, for the smaller-and larger-amplitude oscillations (see Appendix A). If the oscillation amplitude is small as in figure 11, the bubble dynamics is rather insensitive to the value of μ. On the contrary, as the oscillation amplitude gets larger as in figure 12, viscosity will play a more important role; we can clearly observe that the oscillation amplitude is reduced by viscous damping. When it comes to exploring a better fitting of the individual bubble dynamics simulation to the experiment particularly for the case of large-amplitude oscillations, we need to look, instead, for the optimum values of μ and G for each experimental run with a specific bubble nucleus of R 0 .
Finally, we make a comment on the non-spherical oscillations of mode 1 in figure 7(a). According to Murakami et al. (2020), the most unstable mode n of non-spherical bubble oscillations may be predicted based on the nth-order natural frequency (of linearized oscillations), (4.1) This instability occurs when ω = 2πf = 2ω n . In the case in figure 7(a), the equilibrium radius is R 0 = 107.8 μm and the ultrasound frequency is fixed at f = 28 kHz, which enables one to find the most unstable mode n (given surface tension σ , density ρ 0 and shear modulus G of the gel). If using G = 4.0 kPa (i.e. the fitted value from the comparison between the linearized theory and the experiment of the small-amplitude oscillations in figure 5), the most unstable mode is predicted at n = 2, which is inconsistent with the observed result (n = 1). However, following the observation in figure 8 that the nonlinear Rayleigh-Plesset calculation with the doubled shear modulus can yield a better prediction on the resonant bubble radius for the finite-amplitude oscillation, we may use G = 8.0 kPa instead, leading to the most unstable mode at n = 1. This suggests that the linear theory represented by (4.1) may be capable of predicting the most unstable mode appearing in finite-amplitude bubble oscillations, as in figure 7(a).

Conclusions
We studied, with experiments and simulation, nonlinear oscillations of a spherical bubble in a 6 wt% gelatin gel under 28 kHz ultrasound irradiation. Ultrasound intensity and the bubble radius at mechanical equilibrium were treated as parameters in our experiments. The shear viscosity and shear modulus of the gel were obtained by fitting the resonance curve of small-amplitude bubble oscillations as a function of the equilibrium radius to the linearized solutions of the Rayleigh-Plesset model that accounts for the gel's elasticity. We also obtained the resonance curve for larger-amplitude bubble oscillations experimentally and numerically, showing the so-called spring softening effect. To reproduce the experimentally determined resonance curve for larger-amplitude oscillations, it is of great importance to properly select coefficients of the gel viscoelasticity in the Rayleigh-Plesset model, for the resonant radius and peak amplitude are sensitive to the shear modulus and viscosity, respectively. On the other hand, when it comes to exploring a better fit of individual bubble dynamics simulation to the observation of a given experimental run, we need to look for the optimum coefficients of the gel viscoelasticity, instead of those determined from the fitting to the resonance curve for the smaller-amplitude oscillations, as inferred from the supplemental simulations. This implies that the gel structure in the vicinity of the bubble nucleus is expected to be damaged by finite-amplitude oscillations of the bubble.
In our experiments, we observed non-spherical bubble oscillations of mode 1, while we put a focus on spherical bubble dynamics in this study. It is expected that higher-mode shape oscillations appear as gel concentration or ultrasound frequency is reduced (Hamaguchi & Ando 2015;Murakami et al. 2020). Our results showed that the amplitude of volume oscillations of the bubble with the shape oscillations is overestimated by the Rayleigh-Plesset-type simulation, which is a critical drawback when it comes to predicting bubble collapse pressure in the context of therapeutic ultrasound. One future work of interest is to examine effects of shape oscillations on the pressure of bubbles collapsing in viscoelastic solids. in figure 9. We also perform the same simulations based on the linear and nonlinear elasticity models but with the acoustic pressure p A set at 10 times larger than that in figure 7, and compare them in figure 10. These supplemental figures support discussion in § 4.2. Second, we perform the bubble-dynamic simulations as in figure 6(b) (for the case of smaller-amplitude oscillation under f /f R = 1.06) and figure 7(b) (for the case of larger-amplitude oscillation under f /f R = 0.93) with rigidity fixed at G = 4.0 kPa, but with varying viscosity μ; see figures 11 and 12, respectively, for the smaller-and larger-amplitude oscillations. These parametric studies support the discussion in § 4.3.