Parsimonious inertial cavitation rheometry via bubble collapse time

Abstract

The rapid and accurate characterization of soft, viscoelastic materials at high strain rates is of interest in biological and engineering applications. Examples include assessing the extent of tissue ablation during histotripsy procedures and developing injury criteria for the mitigation of blast injuries. The inertial microcavitation rheometry technique (IMR, J. B. Estrada, C. Barajas, D. L. Henann, E. Johnsen and C. Franck, J. Mech. Phys. Solids, 2018, 112, 291–317) allows for the characterization of local viscoelastic properties at strain rates up to 10⁸ s⁻¹. However, IMR typically relies on bright-field videography of a sufficiently translucent sample at ≥1 million frames per second and a simulation-dependent fit optimization process that can require hours of post-processing. Here, we present an improved IMR-style technique, called parsimonious inertial microcavitation rheometry (pIMR), that parsimoniously characterizes surrounding viscoelastic materials. The pIMR approach uses experimental advancements to estimate the time to first collapse of the laser-induced cavity within approximately 20 ns and a theoretical energy balance analysis that yields an approximate collapse time based on the material viscoelasticity parameters. The pIMR method closely matches the accuracy of the original IMR procedure while decreasing the computational cost from hours to seconds and potentially reducing reliance on ultra-high-speed videography. This technique can enable nearly real-time characterization of soft, viscoelastic hydrogels and biological materials with a numerical criterion assessing the correct choice of model. We illustrate the efficacy of the technique on batches of tens of experiments for both soft hydrogels and fluids.

---

1 Introduction

The characterization of viscoelastic soft materials undergoing fast, finite deformations is necessary for a wide range of applications. These include, but are not limited to, the prediction of biological tissue damage due to blunt impact and blast events,^1,2 the design of acoustically-responsive scaffolds for drug delivery,³ and the modeling of non-invasive laser-⁴ and ultrasound-based surgical procedures.^5–7 Notably, the United States Food and Drug Administration (FDA) recently approved the clinical treatment of liver cancer with histotripsy, a novel technique that ablates diseased tissues with ultrasound-induced cavitation.^8,9 A pressing need in histotripsy is an on-the-fly assessment of the degree of therapy completion, which can be reflected through the mechanical response of the treated biomaterial.^10,11 However, soft materials such as hydrogels are challenging to characterize due to their low elastic shear modulus, which ranges from 100 Pa to 1 MPa, and the difficulties of gripping and manipulating the specimens during experiments. To characterize materials with high compliance and, correspondingly, slow shear wave speed, traditional high-strain-rate experiments, such as the Kolsky bar, must be supplemented with pulse shaping, weak signal sensing, and/or other complicating techniques.¹² Furthermore, soft biological tissues often exhibit spatial heterogeneity, increasing the difficulty of measuring their material property distribution with conventional methods that only provide a macroscale average modulus. Bio-inspired material systems fabricated to reproduce these functional gradients are similarly difficult to characterize.

The aforementioned challenges necessitate a technique to locally assess the high strain rate and finite deformation behavior of soft materials. Crosby et al. first developed needle-induced cavitation rheology as an approach to probe the local elastic properties of soft materials.^13,14 A cavity of air or immiscible liquid is injected into the characterized media. The elastic modulus is determined from the pressure and bubble radius at the onset of mechanical instability. This quasi-static approach has been extended in recent years to a ballistic strain-rate regime of approximately 10⁴ s⁻¹ by Milner and Hutchens.^15,16 Cohen and co-workers introduced the capability to cyclically expand and relax the needle-induced cavity at controlled stretch rates,^17,18 enabling finite deformation characterization of viscoelastic materials. The inertial microcavitation rheometry (IMR) technique, introduced by Estrada et al.¹⁹ and improved by others recently,^20–22 accesses a higher range of strain rates by using laser-induced cavitation (LIC) in soft, hydrated materials (i.e., with shear moduli below ∼1 MPa). An ultra-high-speed camera images the bubble kinematics and the viscoelastic properties of the cavitated media are inversely characterized according to an inertial cavitation bubble model^23,24 with refinements accounting for a two-component mixture of bubble content with heat diffusion and mass transfer^25–29 and stress field in the surrounding media.^30–33

IMR inversely characterizes viscoelasticity at strain rates reaching 10³–10⁸ s⁻¹ but has only been successfully applied to characterize nearly transparent materials using ultra-high-speed imaging (rates above 270 000 frames per second). The reliance on the full dynamics acquired from high-speed, bright-field videography of the cavity restricts assessment to experimental systems that produce accurate bubble images, which itself is a product of, e.g., camera sensitivity and the material turbidity. Increasing the exposure time to combat low light throughput works against the maximum frame rate, thus suggesting a need for us to reduce the reliance on the transient dynamics for the characterization of challenging optical systems.

Furthermore, the computational cost of the forward simulation, optimization, and best-fit procedure is restrictive, particularly in a potential desired end-case-usage for near on-the-fly characterization, for example, to assess the extent of diseased tissue removal during histotripsy. Each forward simulation requires about ten seconds. Batch-fitting multiple experiments simultaneously and increasing the number of model parameters cause an exponential increase in the required forward simulations. Hence, we seek to construct an approximate theoretical model that characterizes materials based on just the most essential data drawn from multiple experiments, i.e., maximum radius, quasi-equilibrium radius, and time to first collapse. The potential benefit here of a computationally inexpensive approximate model is twofold. In cases where an approximate characterization of material viscoelasticity is sufficient for predictiveness, the procedure herein represents a rapid method substituting the accurate yet time-consuming IMR procedure. If accuracy remains critical, this procedure complements IMR by vastly paring down the computational space prior to an accurate bubble-dynamics-based inverse characterization.

This style of approximate model can also be extended to applications in which the time to collapse is used to quantitatively describe some system behavior or parameters of interest. For example, the collapse time measured for LIC in the vicinity of agarose hydrogels was compared against the Rayleigh collapse time (a simplified metric assuming the bubble is just a void) by Sieber et al.³⁴ to examine the effect of an elastic boundary. Marsh et al.³⁵ and Ohl et al.³⁶ conducted shock-induced cavitation experiments in water and cervix cell assays, respectively, and approximated average velocity and pressure in the resulting jet flow using the Rayleigh collapse time. These types of analyses could thus be enhanced by our approximate collapse time model accounting for material behavior and other bubble physics.

Herein, we use the modified Rayleigh collapse time approach to develop a strategy for the parsimonious characterization of viscoelastic materials that can be described with up to three-parameter models. In contrast to prior work,¹⁹ this approach enables the use of data from multiple experiments to arrive at a batch-fit solution. The strategy leverages high-fidelity measurements of the maximum bubble radius, the long-term equilibrium bubble radius, and the time from maximum expansion to first bubble collapse. These quantities of interest are distinctly related through the ultra-high-rate elastic and viscous behaviors of soft materials. In Section 2, we present an LIC experiment setup capable of quantifying the time of collapse to an accuracy of approximately 20 ns. The experiments are complemented with an energy balance analysis that approximately quantifies the effects of material viscoelasticity and secondary factors (viz., surface tension, bubble pressure, and dilatational wave speed) on the time to the first bubble collapse. We then introduce the parsimonious inertial microcavitation rheometry (pIMR) procedure enabled by these experimental and theoretical advancements. The consistency of the procedure is verified in Section 3 with synthetic experiments. We demonstrate in Section 4 high-fidelity viscoelastic model parameterization from tens of experiments in viscoelastic liquids and hydrogels, with computational post-processing that takes only seconds. In Section 5, we discuss the implications of the results obtained and the limitations of the proposed strategy. We provide concluding remarks in Section 6.

2 Theory and methods

2.1 Bubble dynamics model

We summarize herein the bubble dynamics model serving as the theoretical basis of the original IMR method. A more thorough discussion of the theory, including its underlying assumptions and regimes of applicability, can be found in Estrada et al.¹⁹ The IMR framework is founded on the classical Keller–Miksis model of bubble dynamics^23,24 and it has been extensively validated for inertial cavitation in nearly-incompressible, viscoelastic, soft materials ranging from polyacrylamide,^19,20 agarose,³⁷ and gelatin^3,22 hydrogels to healthy and diseased human liver tissues.³⁸ Experimental quantification of the deformation field in inertially cavitating hydrogels also confirmed the validity of the bubble dynamics model when the bubble wall velocity is lower than 0.08 times the longitudinal wave speed in the surrounding material.²⁰

Briefly, the bubble dynamics model considers a spherical bubble in an infinite surrounding material environment subjected to a pressure change that causes rapid radial motion, as depicted in Fig. 1. The material outside the cavity is viscoelastic and approximated as nearly incompressible. At the bubble wall, the traction due to the material stress is balanced by the internal pressure of the bubble content and the surface tension. We denote the equilibrium, stress-free radius of the spherical bubble as R₀ and the referential radial coordinate for a material point in the surrounding viscoelastic medium to be r₀ ∈ [R₀, ∞), measured from the center of the bubble to the infinite far field. Due to the balance of mass, the deformed radial coordinate r and velocity v of a material point r₀ at time t in an incompressible medium are

where R(t) is the evolving radius of the bubble. The balance of linear momentum in the radial direction requires thatwhere ρ is the material density of the surrounding material, p the hydrostatic pressure in the material, s_rr and s_θθ the normal radial and normal circumferential components of s, the deviatoric Cauchy stress in the material. There need not be a d/dt term for the density in (2) because we have used continuity of mass to simplify. A perturbation analysis bridging the near- and far-fields of the bubble^23,24,39 leads to a correction of eqn (1) accounting for a finite pressure wave speed c in the material and the energy transfer via outward radial acoustic emission. Then, integrating eqn (2) over r from r = R to r → ∞ and considering the traction boundary condition at the bubble wall results in what is known as the Keller–Miksis equation describing bubble dynamics,where overdots denote derivatives with respect to time t, p_b the internal bubble pressure, p_∞ the far-field pressure, γ the bubble wall surface tension, and S the stress integral defined as

Fig. 1 Schematic representation of the spherical bubble considered in the bubble dynamics and approximate collapse time models. The nearly incompressible, viscoelastic material surrounding the bubble is modeled as a finite deformation, standard linear solid (SLS) described by a ground-state elastic shear modulus G, a viscous shear modulus μ, and a relaxation time scale τ₁. When τ₁ → 0, the SLS model becomes a Kelvin–Voigt model.

We do not simulate the complex plasma physics contributing to the initial growth of the laser-induced cavity. Instead, following a conventional approach for modeling LIC,^{19,25,28,29,40} we assume that the bubble contents and the surrounding medium reaches thermodynamic equilibrium at maximum bubble expansion and model the bubble dynamics from the instance of maximum bubble expansion. When considering surrounding material with history-dependent viscoelasticity, we estimate the initial condition of the stress integral S according to a simplified model of the bubble growth phase, as detailed below in Section 2.1.2 for the Maxwell model. We assume that the bubble contains a mixture of water vapor and other, non-condensible gas components during the rapid bubble dynamics. Mass and heat transfer of the two-part bubble contents are assumed to obey Fick's and Fourier's laws, resulting in a set of PDEs.^26–28 Following earlier works,^19,20,40 we assume that the surrounding material has a sufficiently large heat capacity and thus remains isothermal at an ambient temperature. Numerical solutions to the Keller–Miksis equation coupled with the bubble content equations are obtained with the ode23tb function in MATLAB (The MathWorks, Inc., Natick, MA).

2.1.1 Non-dimensionalization and solution of bubble dynamics model. We follow existing work to non-dimensionalize the governing equations and clarify the interactions between material parameters,¹⁹ and list non-dimensional parameters in Table 1, where R_max is the maximum radius of the bubble, G₁ is the shear modulus associated with the Maxwell element, and is the characteristic velocity. The non-dimensional Keller–Miksis equation describing the evolution of the non-dimensional bubble radius R* is

Table 1 Dimensionless quantities in the Keller–Miksis equation

Dimensional quantity	Dimensionless quantity	Name
t	t* = tvc/Rmax	Time
R	R* = R/Rmax	Bubble-wall radius
R0		Equilibrium bubble-wall radius
c	c* = c/vc	Material wave speed
pb		Bubble pressure
γ	We = p∞Rmax/(2γ)	Weber number
S	S* = S/p∞	Stress integral
G	Ca = p∞/G	Cauchy number
μ	Re = ρvcRmax/μ	Reynolds number
τ1 = μ/G1	De = μvc/(G1Rmax)	Deborah number

---

Unless stated otherwise, we assume ρ = 998.2 kg m⁻³, p_∞ = 101.3 kPa, c = 1484 m s⁻¹, and γ = 0.072 N m⁻¹. For the viscous fluids characterized in Section 4.1, we assume using the rule of mixtures that ρ = 1154 kg m⁻³ for the mixture of water and glycerin and ρ = 1100 kg m⁻³ for the mixture of water and polyethylene glycol (PEG 8000). Assuming a constant temperature of 298.15 K in the surrounding material, the material parameters related to the heat and mass transfer of bubble contents are defined according to Estrada et al.¹⁹

2.1.2 Stress integral in the surrounding medium. The stress integral S* for the viscoelastic constitutive models considered in this work is tabulated in Table 2. We consider finite viscoelasticity constitutive models with stress responses that are additively decomposed into those of three elementary components: a neo-Hookean hyperelastic contribution, a Newtonian viscous contribution, and a Maxwell fading memory viscoelastic contribution. The stress integral for the Kelvin–Voigt viscoelastic models follow our previous work,¹⁹ in which a neo-Hookean hyperelastic spring is arranged in parallel with a Newtonian viscous dashpot. The standard linear solid (SLS) model (sometimes referred to as the Zener model) consists of a neo-Hookean hyperelastic spring parallel to a Maxwell branch.

Table 2 Summary of material stress integrals

Material model	Stress integral relationship S*
Neo-Hookean
Newtonian
Kelvin–Voigt
Maxwell
Standard linear solid (SLS)

---

Assuming that the characteristic time scale of the bubble oscillation is longer than the time scale of the exponential relaxation of a Maxwell material, its stress integral satisfies

Due to the fading memory of the Maxwell material, a non-zero stress integral remains at the end of the bubble growth phase, contributing to the ensuing bubble collapse. This quantity is numerically evaluated by advancing the ODE (6) from the beginning of the growth phase, with initial conditions R = R₀ and Ṙ = Ṙ_i > 0, and Ṙ decreasing to 0 at the end of the growth phase. The fminsearch function in MATLAB is used to iteratively solve for Ṙ_i to minimize |R − R_max| at the end of the growth phase. The value of S at the end of the growth phase is then determined. Heat and mass transfer are neglected in these simulations of the bubble growth phase.

These models, or their modified hyperelastic equivalents, have successfully characterized hydrogels with IMR.^{3,19,20,37,41} In this work, we are primarily interested in the contribution of material viscoelasticity, S*, on the collapse time. We note that the primary non-dimensional parameters of calibration interest, therefore, are the Cauchy number (ground-state elasticity), Ca, the Reynolds number (ground-state viscosity), Re, and/or the Deborah number (relaxation time), De, all defined in Table 1. Thus, in the following sections we distinguish and separately quantify the collapse-time effects from these three material parameters for characterization from those arising from other bubble physics.

2.2 Energy balance analysis and analytical estimates of collapse time

We modify Lord Rayleigh's original analysis to obtain a more accurate prediction of a bubble collapse within hydrogel-like materials. Lord Rayleigh utilized an energy balance approach⁴² with the following four assumptions: (i) the bubble has no contents, (ii) there is no surface tension between the void and the surrounding material, and the surrounding material is (iii) incompressible and (iv) inviscid. Thus, the potential and kinetic energy of the surrounding material dictate the evolution of the bubble radius. Under these conditions, the Keller–Miksis equation (5) simplifies to,where is the non-dimensional liquid pressure . The potential energy of the inviscid liquid surrounding the bubble is the volume integral of the non-dimensional liquid pressure,where is the volume of the bubble. The kinetic energy of the liquid iswhere ρ* is the non-dimensional liquid density (ρ* = 1). The void is assumed to begin at rest, corresponding to an initial kinetic energy of zero. The energy balance is then (4π/3)(R*³ − 1) + 2πR*³Ṙ*² = 0. Isolating the bubble wall velocity as a function of the radius, Ṙ*(R*(t*)), and integrate to the closure of the bubble, we obtain the Rayleigh collapse timewhere Γ[·] is the gamma function. To obtain the dimensional form, we multiply by the characteristic timescale, .

2.2.1 General approach for modified Rayleigh collapse time. We can generalize a Rayleigh-type model as the following equationwhere f* is a sum of different physical phenomena (see Table 3). f* can also be interpreted as a force or a resistance to force. Following the work of Yang et al.,⁴³ for a constant f* an analytical solution for the modified Rayleigh collapse time can be obtained. Thus, we define a time-averaged f* acting on the bubble from the surroundings as

Table 3 Physical phenomena and corresponding functional changes to the Rayleigh–Plesset equation. The right-most column sums to the overall function that transforms eqn (11) into (5) (under the polytropic gas assumption)

Phenomenon	Function modifying Rayleigh–Plesset equation f*
Bubble pressure
Weak compressibility
Surface tension
Material response
Compressibility affecting bubble pressure
Compressibility affecting material response

---

Thus, an ansatz for the change in the liquid potential energy, and corresponding energy balance are,

respectively. Following the procedure of Lord Rayleigh, the approximate bubble wall velocity iswith the general approximate modified collapse time

Since time-averaging is a linear operator, we can write the total collapse time modification to be equal to the following, , where α is indexing different physical effects. We consider the constitutive terms in eqn (5) individually. Since inertia dominates the collapse, the interaction of compressibility with other physical phenomena are second-order and are neglected. Additionally, for simplicity, this analysis will neglect heat and mass transfer in and outside the bubble. Thus, the vapor pressure in the bubble is constant.

The time averaging of f* is similar to linearization in traditional perturbation methods. That is, if any of the forcing terms approach unity, the approximation will break down and exhibit large errors when compared to the exact solution. Neo-Hookean elasticity is an exception, as the leading order elastic forcing term is constant. This exception permits reasonable predictions of the elastic contribution to the collapse time, even for small Ca.

2.2.2 Bubble pressure effects. We assume that there are two primary gases inside the bubble: (1) water vapor and (2) a non-water gas phase. The latter consists of air and vaporized material that diffuses back into the material over time scales much longer than that of inertial collapse. We consider the bubble pressure as the sum of partial pressures of the gases present:⁴⁴ , where κ is the ratio of the heat capacity at constant pressure, C_P, to the heat capacity at constant volume, C_V. Additionally, the water vapor pressure is , and we assume the non-condensible gas to be polytropic, where is the equilibrium bubble pressure.

In the limit R* → 0 (i.e., infinite bubble pressure), the evaluation of the mean value of the bubble pressure is non-convergent. Furthermore, there is no expression available for the minimum radius such that we could obtain a finite integrated result. Since , then and a proportionality constant results through integration of R such that

where [scr B, script letter B] is obtained by numerically solving the exact collapse time integral. The exact collapse time is found by considering the resistive force of the bubble contents preventing collapse. The bubble internal energy, or the reduction of liquid potential energy in the presence of bubble contents, is

For the special (isothermal) case of κ = 1, the bubble internal energy:

Including the bubble internal energy in the energy balance with the liquid potential (8) and kinetic (9) energies, the non-dimensional collapse time is

where

We evaluate eqn (19) for a small (0.01) non-dimensional equilibrium radius. The integral is evaluated by setting the minimum radius to zero. Only the real part of the result is considered, which is equivalent to evaluating the integral from 1 to the minimum radius. For the special cases of κ = 1.4 (isentropic) and κ = 1, we obtain [scr B, script letter B] = 2.1844 and 1.4942, respectively. To be consistent with the initial bubble pressure in IMR,¹⁹ we selected κ = 1 for the pIMR solver.

2.2.3 Weak compressibility effects. In Table 3, the third term of dominates the weak compressibility effect during the primary bubble collapse. The third term is associated with the liquid potential which dominates for most of the collapse phase. Hence, we neglect the other two weak compressibility terms related to the kinetic energy. Thus, , where is the characteristic Mach number. Time averaging for the duration of the collapse and solving explicitly,

2.2.4 Surface tension effects. The surface tension of the water-containing material plays a non-negligible role during the collapse. Time-averaging (see Table 3) we obtain

2.2.5 Neo-Hookean elasticity. Approaching bubble closure, , the neo-Hookean stress integral converges to a constant. Therefore, the modification function and collapse time modification are equivalent for a highly inertial collapse, . Substituting the neo-Hookean expression for into eqn (16) results in the modified collapse time consistent with the result in Yang et al.,⁴³ i.e.,

However, for finite and Ca of [scr O, script letter O](1), the term that is linear in can no longer be assumed to be small. For LIC, and for the majority of the collapse phase, , thus we may neglect the quartic term in in Table 2 and . Accounting for the finite equilibrium bubble radius, the time-averaged f* is

2.2.6 Viscous Newtonian fluid. Unlike the elastic case, the mean value of the modification function for a viscous Newtonian fluid does not converge. For a void to reach closure in a viscous fluid, the strain rate and, therefore, the viscous dissipation becomes infinite. Similar to the approximation in Section 2.2.2, a proportionality coefficient for is to account for an intractable non-zero minimum radius when evaluating f̄*, i.e.,

Here, [capital script C] = [capital script C](Re) to obtain accurate results at smaller Reynolds numbers which are experimentally relevant. Solving eqn (24) for yields

[capital script C] is solved for by balancing the liquid potential (8) and kinetic (9) energies with the viscous dissipation, i.e.,

The bubble wall velocity and exact collapse time are then

respectively. Substituting eqn (15), (16), and (25) into (28), the nested integral on the right hand side is evaluated to obtain an implicit relationship for [capital script C](Re),where ₂F₁(·) is an ordinary hypergeometric function. The right-hand side can be numerically integrated to find the value of [capital script C] for a given Re. For a fast and simple calculation, we approximate the implicit function with a perturbation series where the small parameter is 1/Re such that [capital script C](Re) ≈ [capital script C]₀ + [capital script C]₁/Re + [capital script C]₂/Re². Constants [capital script C]₀, [capital script C]₁, and [capital script C]₂ approximate the implicit function in eqn (29) and are found by numerically integrating and iterating for three separate Reynolds numbers. The Reynolds numbers used for this fitting are 18, 25 and 500. Below a Reynolds number of 18, the numerical integration produces imaginary solutions. The following values approximate eqn (29): [capital script C]₀ = 0.46379, [capital script C]₁ = 0.56391, and [capital script C]₂ = 5.74916.

2.2.7 Kelvin–Voigt viscoelasticity. The Kelvin–Voigt stress integral average is the sum of the viscous and elastic contributions, i.e.,

2.2.8 Maxwell viscoelasticity. The stress integral for the Maxwell model has no closed-form relationship. We approximate the right-hand side of the stress integral ODE in Table 2 to be and obtain the approximate relationship for the stress integralwhere is the unrelaxed stress in the surrounding material at the maximum radius due to the expansion. For De ≪ 1, ; otherwise, the initial stress can alter the subsequent bubble dynamics. The approximation of is described in Section 2.2.10. Evaluating the time-average integral yieldswhere is obtained by substituting eqn (16) into the previous expression. However, the implicit relationship has no analytical solution for . In equation eqn (32), is the actual collapse time that depends on . For simplicity, we approximate the remaining collapse time dependence, , by the Rayleigh collapse time,

2.2.9 Standard linear solid using neo-Hookean elasticity. We consider a material model, the modified standard linear solid or the Zener model,^31,45 comprising a Maxwell element parallel to a neo-Hookean elastic element to be able to describe more complicated finite-deformation viscoelastic material behavior. Since the deviatoric Cauchy stress tensor is a sum of contributions, the stress integral and its time derivative are and , respectively. Similarly, the collapse time modification function is the sum .

2.2.10 Initial stress due to unrelaxed Maxwell element. The energy balance approach is used to approximate the initial, unrelaxed Maxwell stress at the maximum bubble radius. We approximate the growth process as the inverse of the collapse; therefore, the modification functions f̄* switch sign. Additionally, for bubble growth, it is assumed that the fluid is initially stress-free, . Therefore, by using the negative value of , the growth time, and , we approximate the initial Maxwell stress using eqn (31).

To find the growth time, we consider the bubble to be nucleated in a stress-free material at the equilibrium radius with a positive bubble wall velocity such that the correct maximum stretch ratio R_max/R₀ is reached. The non-dimensional energy balance is

where is the unknown initial bubble wall velocity and the average of the Rayleigh–Plesset modification function during the growth phase. The initial bubble wall velocity is obtained by setting the non-dimensional current bubble radius and bubble wall velocity to 1 and 0, respectively. Solving for the initial bubble wall velocity yields:

The current bubble wall velocity is obtained by substituting eqn (35) into eqn (34) and taking the positive root for bubble growth:

The growth time is then

If , then the growth time approximation has an imaginary, unphysical contribution due to averaging the forcing during the growth phase. The two physical effects that can produce these imaginary solutions are elasticity and surface tension. Neo-Hookean elasticity and surface tension produce unphysical results for Ca < 5/2 (G > 50 kPa) and We of the same order as Ca, respectively. The shear moduli observed in this work were all below 25 kPa. For We of [scr O, script letter O](1), the maximum radius would be much smaller than is experimentally relevant in this work. Small Deborah numbers (<0.1) result in the initial stress prediction deviating from the iterative method result described in Section 3 with relative errors above 50% (data not shown). However, the accuracy of the initial stress in this regime is inconsequential, as the corresponding collapse time modification factor is very small. For the collapse time prediction, parameter values yielding a f̄* > 1 would lead to weak oscillations not inertial bubble collapse.

2.3 Experimental methods

The laser microcavitation experiments follow the general LIC procedure of Estrada et al.¹⁹ with two main advancements of (i) shadowgraph and ghost imaging and (ii) incident beam shaping⁴⁶ (see Fig. 2). These are described in detail in Appendix A and summarized here.

Fig. 2 The experimental setup to generate, record, and profile single laser-induced microcavitation (LIC) bubble events in soft materials. The setup uses a combination of a class-4, frequency-doubled Q-switched 532 nm Nd:YAG pulsed laser, a high-speed imaging camera, and a spatial light modulator. The time of the first bubble collapse is estimated according to the shock wave, which was visualized by shadowgraph and ghost imaging techniques.

Single LIC bubble events are generated in soft materials using a pulsed, frequency-doubled (532 nm), Q-switched Nd:YAG laser. The pulse energy is user-defined and was on the order of 1–10 mJ for the experiments. A diffraction-limited focusing objective condenses the laser pulse into a beam waist to approximately 4 μm in diameter. A second objective is oriented orthogonal to the imaging plane for the purpose of verifying bubble sphericity. A spatial light modulator is used to tune higher-order beam asymmetry to create spherical bubble events. We record the microcavitation event at 1 million frames per second (Mfps) using a Shimadzu HPV-X2 (Tokyo, Japan) ultra-high-speed imaging camera. A backlight laser strobe fires synchronously with the camera and is sent to the bubble event as parallel light. This light enables shadowgraph imaging, a mode related to Schlieren imaging that permits the visualization and measurement of emitted shockwaves. We strobe the shadowgraph backlight twice per frame, improving our estimate of the collapse time using the shock speed (found by locating two shocks on one frame) and the minimum radius estimate.

Water–glycerol solutions were prepared by combining glycerol (G33-1, Thermo Fisher) with DI water at a volumetric ratio of 60% glycerol to 40% DI water. The solution was subsequently placed on a magnetic stirrer (Isotemp SP88854200, Thermo Fisher) and stirred continuously at room temperature for a duration of 30 minutes. Water–PEG mixture samples were produced by mixing polyethylene glycol 50% (w/v) of molecular weight 8000 (PEG 8000; Avantor 101443-878) with DI water in a proportion of 80% PEG by volume. The blended mixtures were poured into glass-bottomed, 35 mm diameter Petri dishes up to roughly 2 mm of depth. These prepared samples retained a liquid state with no signs of heterogeneity. The low-frequency shear moduli of the water–PEG mixture samples were measured using a TA Instrument ARES-G2 rotational rheometer (New Castle, DE) equipped with a 40 mm diameter stainless steel 2°-angle cone plate fixture and a flat base. Dynamic loading was applied at a frequency of 1 rad s⁻¹ with the maximum strain amplitude of 0.04 rad.

Polyacrylamide (PAAm) gels for characterization purposes were prepared at concentrations of 5%/0.03% and 10%/0.06% acrylamide/bisacrylamide (v/v) according to previously developed protocols.^19,47 The PAAm gels were cast in square 5 mL polystyrene spectrophotometer cuvettes and cured for 45 min prior to characterization.

2.4 Parsimonious inverse characterization based on bubble collapse time

Past studies using IMR have found that adjusting the laser energy can modulate the maximum radius of the bubble, while the amplification factor of the initial bubble expansion, Λ_max = R_max/R₀, is weakly sensitive to laser energy.⁴⁸ Thus, we perform LIC experiments at various laser energy levels on a material and tune the parameters appearing in eqn (5). Our experiments traverse the {R_max, Λ_max} space for a constant set of dimensional viscoelastic parameters and collect . As illustrated in Fig. 3 for a Kelvin–Voigt model example, reflects the combined effect of the collapse modification functions reviewed in Section 2.2, which varies with R_max and Λ_max.

Fig. 3 Combined effect of viscoelasticity, bubble content pressure, weak compressibility, and surface tension on the bubble collapse time in a Kelvin–Voigt material with {G = 10 kPa, μ = 0.1 Pa s} across typical range of R_max and Λ_max in LIC experiments.

We solve for the viscoelastic model parameters that minimize the difference between the collapse times approximated by the energy balance analysis, t^Approx_c, and those that were experimentally measured, t^Expt_c. We refer to this inverse characterization method as the parsimonious inertial microcavitation rheometry technique (pIMR).

Specifically, we use a cost function,

that quantifies the agreement between the collapse time measured experimentally and the approximated value for a set of trial parameters {G, μ, τ₁} according to the energy balance analysis using n number of experiments. The cost function can be interpreted as the order of mean square relative error between the measured and predicted collapse time. Using the fminsearch function in MATLAB, which implements a Nelder–Mead direct search process,^49,50 an optimal set of viscoelastic parameters is then determined to minimize ψ.

To analyze the precision of the inverse characterization solution, we also introduce the normalized cost function for experimental data, [small psi, Greek, circumflex] = ψ − ψ₀, where ψ₀ is the minimized cost function for a given type of viscoelastic model, corresponding to the optimal solution found by the direct search algorithm. The normalized cost [small psi, Greek, circumflex] is equal to zero at the optimal solution, whereas the positive-valued [small psi, Greek, circumflex] elsewhere reflects how far the solution is from being optimal.

Furthermore, we seek to encourage the parsimony of the inversely-calibrated viscoelastic model type and minimize the number of parameters used. This could be achieved, for example, through a F-test-based criterion that discourages the addition of a model parameter that does not lead to a large enough decrement of ψ₀. In practice, a user could decide to penalize a model multiplicatively based on the added number of terms^51–53 or use a least absolute shrinkage and selection operator (LASSO) regression.⁵⁴ In the present work, we simply report cost function decrement Δψ₀ to reflect the parsimony of the constitutive model. For the specific models we consider, we report the cost function decrements

where ψ_0,NH, ψ_0,Nt, ψ_0,KV, ψ_0,SLS are the minimized ψ corresponding to the neo-Hookean, Newtonian, Kelvin–Voigt, and SLS models, respectively. If a model type results in a decrement below a threshold value of Δψ₀, we can consider it to be over-fitting. In our analysis presented below, a threshold value of 0.5 is considered for an illustrative purpose. A user may modify the choice of threshold depending on the type of material characterized and the relative amount of noise in the experiment data.

Parsimonious inertial microcavitation rheometry procedure.

Data generation: Follow methods in Section 2.3 to obtain R(t) and t_c for n experiments.

From dataset, R_max, R₀, calculate and then given material properties: c, γ, ρ then

 Calculate M_c and We from dataset of R_max

 Calculate and

Provide initial guess for constitutive material properties, μ₀, G₀, and τ_1,0

For each constitutive model, calculate an initial loss, ψ_i while ψ ≠ ψ₀ do

 Evaluate total f̄* and calculate t^Approx_c for each model

 Compute updated loss, via (38)

  Update parameters via Nelder–Mead step (fminsearch in MATLAB)

end while

Compute both Δψ₀ to determine model via (39)

*Optional

 To verify model choice, compute RMSE values between R^Expt and R^Sim for all models

 using μ, G, and τ₁ from pIMR

 if refine material constitutive parameters further do

  Traditional IMR fitting procedure via least squares minimization on optimal model

  from verification, see ref. 19

 end if

Result: Estimate material model and associated properties, μ, G, and τ₁

3 Consistency check of pIMR

To check the consistency of the pIMR procedure, we use it to recover input viscoelastic models from synthetic experiments. Using the bubble dynamics model and stress integral evaluation procedure described in Section 2.1, we simulate R(t) corresponding to n = 36 pairs of R_max ∈ [100, 400] μm and Λ_max ∈ [5, 9]. The R_max and Λ_max values are generated with the Latin hypercube sampling method assuming a uniform distribution.⁵⁵ The simulated collapse time is then used to inversely calibrate viscoelastic models with pIMR. The results are presented in Table 4.

Table 4 Calibrated viscoelastic parameters, minimized cost function, and cost function decrement from synthetic experiments. (Numerical values below 10⁻⁵ are reported as ∼0)

Material	Model	G (kPa)	μ (Pa s)	τ1 (μs)	ψ0	Δψ0
Synthetic NH, G = 10 kPa	NH	10.44	—	—	−6.09	—
	Newtonian	—	∼0	—	−1.47	—
	KV	10.44	∼0	—	−6.09	0.00
Synthetic Newtonian, μ = 0.1 Pa s	NH	∼0	—	—	−1.83	—
	Newtonian	—	0.096	—	−4.85	—
	KV	∼0	0.096	—	−4.85	0.00
Synthetic KV, G = 10 kPa, μ = 0.1 Pa s	NH	5.07	—	—	−2.79	—
	Newtonian	—	∼0	—	−1.94	—
	KV	10.11	0.095	—	−5.50	2.71
	SLS	10.38	0.116	0.506	−5.78	0.28
Synthetic SLS, G = 10 kPa, μ = 0.1 Pa s, τ = 1 μs	NH	6.56	—	—	−3.41	—
	Newtonian	—	∼0	—	−1.80	—
	KV	9.17	0.052	—	−5.16	1.75
	SLS	10.08	0.096	1.90	−5.97	0.81

---

If additional constitutive model parameters are rejected when Δψ₀ < 0.5, for example, pIMR can correctly identify the type of constitutive model used in the synthetic experiments. Compared to the input values, the elastic and viscous shear moduli, G and μ, are recovered to within an accuracy of 5%. This is well within the confidence interval commonly reported for hydrogels characterized by IMR,^19,20,37,41 confined and unconfined compression,^56–58 and indentation.⁵⁶ For the SLS model, the relaxation time scale τ₁ is recovered to within a factor of two of the input value. This accuracy is acceptable since τ₁ contributes to the viscoelastic stress through an exponential relaxation function.

In Appendix B, additional synthetic experiments with n < 36 and artificially perturbed collapse time data are considered. We find that the accuracy of collapse time measurement in our LIC experiment setup is sufficient to ensure the stable performance of pIMR. Although it is theoretically possible to calibrate a constitutive model with m parameters using data from n = m LIC experiments, this makes the calibration results more susceptible to the inherent discrepancy between the bubble dynamics model and the approximate collapse time model presented in Section 2.2. A larger sample size n is encouraged for the accurate calibration of viscoelastic model type and parameters.

We also verified that the assumption of isothermal bubble content made in Section 2.2.2 has a negligible effect on the inverse calibration results. Using collapse time data from synthetic experiments with isothermal bubble content, without heat and mass transfer, pIMR recovered viscoelastic models matching those shown in Table 4.

4 Inverse characterization of viscoelastic materials

The proposed pIMR procedure is applied to the inverse characterization of viscoelastic materials from LIC experiment results. The calibrated viscoelastic model parameters are summarized in Table 5. The range of R_max and Λ_max spanned in the experiments are shown in Fig. 7(b).

Table 5 Inversely characterized viscoelastic parameters from LIC experiments

Material	Technique	Model	G (kPa)	μ (Pa s)	τ1 (μs)	ψ0	Δψ0
Water–glycerol	pIMR (n = 11)	Newtonian	—	0.012	—	−2.94	—
	pIMR (n = 11)	KV	4.07	0.110	—	−3.03	0.09
	IMR	Newtonian	—	0.066	—	—	—
	Capillary viscometer^59	Newtonian	—	0.010	—	—	—
Water–PEG	pIMR (n = 20)	Newtonian	—	0.223	—	−1.72	—
	pIMR (n = 20)	KV	∼0	0.223	—	−1.72	0.00
	IMR	Newtonian	—	0.151	—	—	—
	Shear-plate rheometry	Newtonian	—	0.122 ± 0.005	—	—	—
PAAm, 5/0.03%	pIMR (n = 52)	NH	3.11	—	—	−3.19	—
	pIMR (n = 52)	KV	6.52	0.109	—	−3.28	0.09
	pIMR (n = 52)	SLS	18.42	0.731	6.24	−3.34	0.06
	IMR	KV	5.01	0.145	—	—	—
	Static compression^19	NH	0.461 ± 0.004	—	—	—	—
PAAm, 10/0.06%	pIMR (n = 39)	NH	10.11	—	—	−3.18	—
	pIMR (n = 39)	KV	14.49	0.130	—	−3.29	0.11
	pIMR (n = 39)	SLS	21.31	0.538	6.03	−3.30	0.01
	IMR	KV	12.02	0.115	—	—	—
	Static compression^19	NH	2.97 ± 0.06	—	—	—	—

---

4.1 Characterization of viscous fluids

We characterized mixtures of water and glycerol, with a 60% v/v glycerol concentration. The mixtures are expected to exhibit viscous fluid behavior with negligible elasticity. A total of 11 LIC experiments were considered, with R_max ranging from 186.8 μm to 345.7 μm and Λ_max ranging from 3.55 to 4.35. The inverse fitting of the Newtonian model resulted in a viscosity of 0.012 Pa s, which is consistent with previously reported values.⁵⁹ Fitting further with the Kelvin–Voigt model resulted in a minimal cost function decrement. Fig. 4(a) shows the approximated collapse time t^Approx_c for the calibrated model versus the measured collapse time t^Expt_c of each experiment. We observe that, on average, t^Expt_c is slightly larger than the predicted value for an inviscid material. This suggests the dominance of material viscosity over elasticity during the bubble collapse. Due to the low value of μ, the effect of the material viscosity is not pronounced in Fig. 4(a). However, the calibrated viscosity of the Newtonian model indeed corresponds to a minimization of ψ and an improved agreement between t^Expt_c and t^Approx_c.

Fig. 4 Comparison of measured vs. predicted collapse time for (a) the 60% (v/v) concentration water–glycerol mixture characterized, with n = 11 samples and (b) the 80% (v/v) concentration water–PEG mixture characterized, with n = 20 samples. The collapse time during LIC is larger than the predicted value for an inviscid material, suggesting the dominance of viscosity over elasticity during the collapse process.

Next, we characterized uncured mixtures of water and PEG 8000, with an 80% v/v PEG concentration. A total of 20 LIC experiments were performed, with R_max ranging from 103.7 μm to 342.9 μm and Λ_max ranging from 4.73 to 7.10. The inverse fitting of the Kelvin–Voigt model converged to a Newtonian model with minimal elasticity. Fig. 4(b) shows the approximated collapse time t^Approx_c for the calibrated model versus the measured collapse time t^Expt_c of each experiment. If the water–PEG material were Newtonian, as in the case of the water–glycerol experiments shown in Fig. 4(a), the predicted and experimentally observed collapse times would fall along the dashed red line corresponding to an agreement between the measurements and predictions. However, the Newtonian pIMR model overpredicts the collapse time for experimental collapse times less than 25 μs, and underpredicts for longer times. Therefore, we hypothesize that the fluid is exhibiting non-Newtonian behavior in the high-strain rate regime.

Using the IMR technique, a Newtonian model with μ = 0.151 Pa s was found to minimize the offset between the normalized bubble history {t*, R*(t*)} recorded experimentally and simulated by the bubble dynamics model, up to the third oscillation peak. Fig. 5(a) shows R(t) for a typical experiment with the simulated bubble dynamics and the inversely characterized constitutive model parameters. As we expect, the optimal parameters of a Newtonian model calculated via pIMR produce dynamics that closely match the bubble collapse time, with an error of 0.29 μs (relative error: 0.94%). The IMR-calibrated model reproduced the post-collapse bubble dynamics more accurately than pIMR but failed to capture the correct collapse time, as shown in Fig. 5(a).

Fig. 5 Bubble dynamics corresponding to representative experiment data (hollow squares) and the inverse characterization solutions: (a) water–PEG 8000 mixture, (b) PAAm gel with 5/0.03% (v/v) acrylamide/bisacrylamide concentration. (c) Decomposition of collapse modification function f̄* for the representative experiments considered in (a) and (b). Relative to the weak compressibility (, teal) and the surface tension (, pink), the bubble pressure (, orange) and the viscoelastic stress (, yellow) have stronger modifying effects on the bubble collapse time.

The least-squares fitting method employed by IMR obtains agreement between simulation and experimental data of the entire transient bubble dynamics. As a result, individual time events, such as the primary bubble collapse time, can be inaccurate. Thus, to accurately reproduce the post-collapse dynamics of non-Newtonian fluids, shear-dependent viscosity models, e.g., Carreau model, are needed for pIMR. While this behavior is not the focus of this work, it does warrant further investigation.

4.2 Characterization of polyacrylamide gels

We characterize PAAm gels with two different concentrations of acrylamide/bisacrylamide. This class of material has been characterized with IMR in past studies^19,40 and exhibited viscoelastic behaviors that were captured well by the Kelvin–Voigt model.

A total of 52 LIC experiments were performed on specimens with an acrylamide/bisacrylamide concentration of 5/0.03% (v/v), with R_max = 218.0–401.3 μm and Λ_max = 6.49–8.46. The history of R(t) for a representative experiment is shown in Fig. 5(b) with the simulated bubble dynamics of the calibrated models. The bubble dynamics of the pIMR neo-Hookean and Kelvin–Voigt solutions matched the collapse time within 0.087 μs (relative error: 0.26%) and 0.063 μs (relative error: 0.19%), respectively. The IMR-calibrated Kelvin–Voigt model overestimates the collapse time by 0.84 μs (relative error: 2.5%). If we were to reject a constitutive model when Δψ₀ < 0.5, for example, the calculated Δψ₀ suggests that the neo-Hookean model suffices to describe the scaling of bubble collapse time. However, the Kelvin–Voigt model clearly reproduces the post-collapse bubble dynamics more accurately in Fig. 5(b).

A total of 39 LIC experiments were performed on gels with an acrylamide/bisacrylamide concentration of 10/0.06% (v/v), with R_max = 215.2–416.3 μm and Λ_max = 5.67–6.76. Again, cost function decrement Δψ₀ suggests that the collapse time scaling is described well by the neo-Hookean model. Consistent with the IMR calibration results, pIMR suggests that the elastic modulus increases with the concentration of acrylamide/bisacrylamide, while the viscosity changes minimally between the two types of specimens.

5 Discussion

The inverse characterization results in Section 4 show pIMR estimating finite deformation viscoelastic model parameters across a batch of experiments with different R_max and Λ_max. For the 52-sample batch of PAAm gel (5/0.03% (v/v) acryl/bis) experiments, the optimal model type and parameters for all samples were determined within 1 second of computation on a workstation (Intel Core i7 14700K). Using the IMR bubble dynamics model, approximately 10 seconds of computational time are required to simulate the bubble dynamics up to the fourth peak of oscillation for each set of input parameters describing the material viscoelasticity and the bubble's initial and equilibrium conditions in each experiment. The computational cost is amplified as the simulation is repeated for combinations of input parameters. In the context of an end-case-usage to characterize soft tissues on-the-fly during medical procedures, the reduction in computational cost from pIMR corresponds to a decreased surgical time and an increased rate of patients treated. While IMR can be accelerated through alternative approaches such as Bayesian optimal design,^60,61 such a procedure still necessitates forward simulations of bubble dynamics. Inverse characterization considering the full bubble dynamics is unable to pare down the parameter space as rapidly as the collapse-time-based, analytical model presented here for pIMR.

The estimation of viscoelastic properties from collapse time also reduces the requirements that IMR previously placed on the optical turbidity of the characterized material. With a decreased frame rate and an increased exposure time per frame, bright-field videography can be used to measure the maximum and equilibrium radii of the bubble in an optically turbid material. Since the bubble collapse coincides with the emission of shock waves, its occurrence can be captured with methods other than the optical strategy introduced in Section 2.3. Integrated circuit piezoelectric transducers are commonly used in shock tube^35,62,63 and Kolsky bar^12,64 experiments to detect pressure spikes during high strain rate deformation of materials. Past studies of laser- and ultrasound-induced cavitation have used hydrophones to acquire acoustic signals and identify the occurrences of shockwave-emitting collapse events.^28,65 Using custom-built histotripsy arrays with receive capable elements, Sukovich et al. have demonstrated the experimental quantification of the time lapse between the nucleation and the first collapse of ultrasound-induced cavitation in ex vivo porcine and bovine tissues.^10,11,66 These acoustic techniques can be feasibly integrated into a LIC experiment setup for the inverse characterization of optically turbid materials.

The inverse characterization of PAAm gels suggest that, at the length scale of LIC experiments, the material elasticity has a stronger contribution to the bubble collapse time than the material viscosity. This agrees with past studies concluding that the first collapse of LIC in hydrogels are dominated by inertial and elastic effects.^20,43 In Fig. 6(a) and (b), the approximated collapse times t^Approx_c for calibrated models are plotted against the measured collapse time t^Expt_c of each experiments in the two types of PAAm gels. For each LIC experiment, t^Expt_c is shorter than what is predicted for an inviscid fluid, indicating that the ground-state elasticity is more dominant than the material viscosity during the bubble collapse. In Fig. 5(c), the collapse modification function f̄* in each representative experiment is decomposed into the contributing parts listed in Table 3. For both the water–PEG mixture and PAAm gel cases, the decomposed values indicate that the bubble pressure and the viscoelastic stress are the main factors modifying the bubble collapse time. The effects of the material compressibility and the surface tension are weaker in comparison. As shown in Table 5, the neo-Hookean model sufficiently decreases ψ₀ and seems to be the optimal choice of constitutive model for the PAAm gel. However, a close inspection of Fig. 6(a) and (b) reveals that the calibrated neo-Hookean model underestimated the collapse time in experiments with small R_max and overestimated the collapse time in experiments with large R_max. The addition of material viscosity improved the agreement between the measured and predicted collapse time values since increases in magnitude with smaller Re, which is inversely proportional to R_max when μ is constant, as discussed in Section 2.2.6.

Fig. 6 Characterization of PAAm gels with pIMR. Comparison of measured vs. predicted collapse time for (a) 5/0.03% (v/v) acrylamide/bisacrylamide and (b) 10/0.06% (v/v) acrylamide/bisacrylamide for inviscid, neo-Hookean, Kelvin–Voigt, and standard linear solid models. Contours of normalized cost function [small psi, Greek, circumflex] corresponding to Kelvin–Voigt parameters {G, μ} for (c) 5/0.03% (v/v) acrylamide/bisacrylamide and (d) 10/0.06% (v/v) acrylamide/bisacrylamide experiments. The opposing effects of elastic modulus G and viscous modulus μ on the bubble collapse time are reflected in the slope of the [small psi, Greek, circumflex] space.

The LIC experiments of PAAm gels surveyed ranges of R_max with a ratio of ∼2 between the upper and lower bounds. The Kelvin–Voigt parameters, G and μ, calibrated with pIMR have relative errors within 30% of the IMR result. However, the normalized cost function spaces shown in Fig. 6(c) and (d) reflect a lower precision in the calibration of μ compared to G. For example, in the case of the 10/0.06% (v/v) acrylamide/bisacrylamide experiments, the region of [small psi, Greek, circumflex] < 0.1 spans 11.1 kPa ≤ G ≤ 18.2 kPa and 0.03 Pa s ≤ μ ≤ 0.22 Pa s, corresponding to upper-to-lower bound ratios of 1.6 and 7.3, respectively. The precision of the calibrated viscosity can be improved by surveying as broad of a range of bubble sizes as is experimentally feasible. Because the viscous forcing function depends on Re, which is linearly proportional to R_max, maximizing the experimental range of R_max would distinguish this effect of material viscosity on the bubble collapse time. In our current setup, practical experimental considerations bound the maximum (in this case, due to the finite cuvette size) and minimum values (e.g., due to the camera frame rate) of R_max. The range of R_max for LIC experiments could, in general, be broadened with longer-focal-length objectives permitting the generation of larger bubbles in larger samples, or sub-nanosecond laser pulses to produce more reliable small bubble events. In the case of more complicated design spaces than R_max alone (such as a variety of external pressures or initial stretch ratios), the pIMR method could be performed in tandem with a recently developed Bayesian optimized experimental design procedure by Chu et al.^60,61 currently only employing bubble dynamics forward simulations. In this scenario, the pIMR approach further stands to speed up characterization by using the analytical model for informing the next best experiment to run for maximum information gain.

In addition to the initial collapse time considered in the present work, other measurable parameters or constitutive models with additional effects (e.g., non-Newtonian behavior) may be harnessed for a more effective parsimonious characterization of viscoelastic models. For example, Fig. 5(b) shows that the neo-Hookean and Kelvin–Voigt solutions from pIMR lead to bubble dynamics that diverge more discernibly from each other after the initial collapse. Similarly, for the case of a water–PEG mixture, shown in Fig. 5(a), the difference between the Newtonian fluid models calibrated via IMR and pIMR becomes clearer in the post-collapse bubble dynamics. Comparing the diverging solutions between pIMR and IMR past the primary bubble collapse, a hybridized approach containing the speed of pIMR and accuracy of IMR is warranted. In such a hybrid setting, pIMR and IMR are complementary to each other. That is, pIMR is used to identify the constitutive model and estimates of the associated material parameters, which may be further refined using an IMR inverse characterization within a pared-down range of constitutive model parameters.

6 Conclusions

We present the pIMR technique, a parsimonious enhancement of the IMR technique that rapidly characterizes the local viscoelastic properties of soft materials from laser-induced cavitation experiments. This new procedure is possible due to experimental advancements in estimating the collapse time of a laser-induced cavity, coupled with a theoretical energy balance analysis. We make an ansatz to a modified potential energy through averaging effects within the Keller–Miksis equation. This ansatz allows the collapse time approximation to include viscoelastic parameters, surface tension, bubble pressure, and finite wave speed. In our approach, we do not introduce empirical fitting parameters in the energy balance analysis to improve its agreement with the bubble dynamics model. These approximate models for the collapse time were shown to perform well in predicting the collapse time from simulations of the Keller–Miksis equation over a parameter space that is experimentally relevant to inertial microcavitation within soft materials.

The proposed procedure successfully pares down the space upon which we seek the global optima of viscoelastic model parameters. Using a cost function ψ that quantifies the agreement between the measured and predicted collapse time, our procedure identifies the simplest type of constitutive model and the optimal values of model parameters. Experimental characterization of viscous fluid and hydrogel specimens resulted in optimized Newtonian and Kelvin–Voigt parameters, respectively, that closely matched the results of the IMR procedure while reducing the computational cost of post-processing from more than an hour to a few seconds.

Our LIC experiments in viscous fluids and soft hydrogels revealed that the dominating mechanisms during the first collapse of the bubble do not necessarily dominate during the ensuing bubble dynamics. For the case of PAAm gels, a neo-Hookean hyperelastic model suffices to reproduce the bubble collapse, while the post-collapse kinematics is strongly influenced by the material viscosity and described better by a Kelvin–Voigt model. We envision that this issue can be addressed via an inverse characterization procedure considering additional observable parameters in the post-collapse bubble dynamics, additional physics in the models (e.g., non-Newtonian behavior), or through coupling pIMR with IMR.

While the present work only considers viscoelastic models with three or fewer parameters in part due to potential non-uniqueness of solutions involving collapse time alone, the procedure of finding the corresponding modification functions can be straightforwardly extended to viscoelastic models with other non-linear elastic and non-Newtonian fluid behaviors. For example, by recording data over a range of stretch ratios as in ref. 48 and incorporating higher-order elastic spring elements,²⁰ the non-linear elastic response can be decoupled from the shear modulus. Quantification of higher-order non-linear behavior is a potentially impactful future use case of this method, as an apparent assessment need in histotripsy is the quantification of the degree of therapy completion via acoustic emissions.^10,11 Haskell et al.¹¹ found an increase in the time from bubble initiation to collapse during the course of therapy, which is qualitatively due to conversion of elastic biomaterial to an ablated viscous liquid. The method presented herein could quantitatively describe the material mechanics during the course of the therapy, and hence, the time to therapy completion. We thus anticipate pIMR to be a useful tool in establishing mechanics-based therapy guidelines for different prospective tissue applications.

Author contributions

ZZ and SR: conceptualization, methodology, software, validation, formal analysis, investigation, writing – original draft, writing – review and editing, visualization. BAA: methodology, software, validation, formal analysis, investigation, writing – original draft, writing – review and editing. DF: methodology, validation, formal analysis, investigation, writing – original draft, writing – review and editing. SHB: writing – original draft, writing – review and editing, supervision, funding acquisition. JY: conceptualization, writing – original draft, writing – review and editing, supervision, project administration, funding acquisition. MR and JBE: conceptualization, methodology, validation, writing – original draft, writing – review and editing, supervision, project administration, funding acquisition.

Conflicts of interest

The authors have no known conflicts of interest associated with this publication, and there has been no significant financial support for this work that could have influenced its outcome.

Data availability

Supplementary data and codes are available at Michigan's Deep Blue Repository. See DOI: https://doi.org/10.7302/fpcw-1173.

The code for pIMR is available at: https://github.com/InertialMicrocavitationRheometry/parsimonious_IMR.

Appendices

A Experimental methods

Our setup generates, records, and profiles pulses of single LIC bubble events in soft materials using a combination of a pulsed, Q-switched, user-adjustable 1–25 mJ, frequency-doubled 532 nm Nd:YAG laser (Continuum Minilite II, San Jose, CA) and a high-speed imaging camera (HPV-X2; Shimadzu, Kyoto, Japan). The setup is triggered using an 8-channel pulse/delay generator (Model 577; Berkeley Nucleonics, San Rafael, CA) according to a customized pre-programmed pulse sequence. The pulse sequence was validated using an oscilloscope (P2025; Berkely Nucleonics). Sequential triggering signals fire two single pulses: the first triggers the laser's flash lamp, and the second fires the Q-switch. The last two triggering signals are sent to a beam profiler (BC106N-VIS; Thorlabs) and the high-speed camera. The backside of the sample is illuminated with the aid of a 640 nm monochromatic ultra-high-speed strobed diode laser (Cavilux Smart UHS; Cavitar, Tampere, Finland). The high-speed camera sync-out signal triggers the illumination laser. The laser beam/pulse was aligned to the back apparatus of a 10×/0.25 high-power microspot focusing objective (LMH-10X-532; Thorlabs, Newton, NJ) using three reflective broadband dielectric mirrors (BB1-E02; Thorlabs, Newton, NJ), three short-pass dichroic mirrors, a beam-sampler, and a spatial light modulator (SLM) (Holoeye, Berlin, Germany). The first dichroic mirror (DMSP605; Thorlabs, Newton, NJ) is used for the beam alignment in conjunction with a continuous exposure collimated laser diode module (CPS635R; Thorlabs, Newton, NJ). A 2× fixed magnification beam-expander (GBE02-A; Thorlabs, Newton, NJ) helps distribute the collimated beam on a larger area and minimizes any potential damage to the SLM and focusing objective lens at the back aperture. The second high-pass dichroic mirror (DMSP550; Thorlabs, Newton, NJ), which has a cutoff wavelength of 550 nm, was used to filter infrared wavelengths and discard them into a beam-block (LB2; Thorlabs, Newton, NJ). The visible beam is then reflected onto a spatial light modulator, allowing for higher control over the last pulse shape and energy. Last, the beam is split before it reaches the focusing objective using a beam sampler lens (BSF10-A; Thorlabs, Newton, NJ). Approximately 0.5% of the split beam is reflected towards a beam profiler (BC106N-VIS; Thorlabs, Newton, NJ) to assess the pulse quality and measure its energy. The remaining 99.5% of the beam continues to the focusing objective through the third dichroic mirror (DMSP550; Thorlabs, Newton, NJ), which also has a cutoff wavelength of 550 nm, allowing the cavitation laser (532 nm) to pass while reflecting the illumination laser light (640 nm). The focusing objective focuses the beam at the microcavitation imaging plane.

The microcavitation event is performed at 1 million frames per second (Mfps) using a Shimadzu HPV-X2 (Tokyo, Japan) high-speed imaging camera, illuminated by CAVILUX Smart UHS (Tampere, Finland) laser, and through both, the cavitation objective and an Olympus Plan 10X-0.25 Achromat imaging objective (RMS10X; Thorlabs, Newton, NJ). The data is analyzed using our in-house Matlab image processing code. To measure the wave speed in the medium, we deployed two imaging techniques simultaneously: laser shadowgraph⁶⁷ and ghost imaging.⁶⁸ Shadowgraph imaging is performed by manipulating the backlighting path to capture density variation due to the compressive shockwave. The physical location of the pressure wave is then estimated during the bubble's cavitation and collapse. Ghost imaging is achieved by triggering the strobed backlight a user-defined number of times per camera exposure, usually 2 or 3 per frame.

B Additional consistency check of pIMR

In this section, we first present a brief summary of how the collapse time models perform when compared to simulations, then we show additional synthetic experiments to verify the consistency of pIMR. We consider a Kelvin–Voigt material with G = 10 kPa and μ = 0.1 Pa s.

We evaluate performance of the collapse time models by computing a relative error between the predicted to the simulated values. We simulate and compare the effects from each non-dimensional parameter individually to determine ranges of the parameters where the theory holds. As such, we numerically solve Rayleigh's equation with only one modification on the right hand side, that is,

where α indexes effects from each non-dimensional parameter. To quantify ranges for the parameters in which the theory is valid, we set a threshold on the relative error in order to determine these ranges. That is, if the relative error between the simulated collapse time and the predicted collapse time is <1%, then the theory holds.

Table 6 shows the results of quantifying regions of validity in the parameter space for the approximations.

Table 6 Non-dimensional parameters with associated theoretically and experimentally valid regimes

Parameters	Theory-valid range	Max theory-valid value of f̄*	Experimental range
a The elastic effect due to a neo-Hookean element does not obey the same linearization rules as the other models. This behavior is because the leading order term in the modification function, f*, is constant. As shown in ref. 43, for a void where , this is the exact modification (f̄*) to the collapse time.b No improvement was observed from including relaxation in the constitutive models when calibrating using IMR. De is of [scr O, script letter O](1) to observe the effect of relaxation.c True for Re > 70. If the Reynolds number is smaller, relative errors >1% are observed for De between 0.03 and 13 (for Re = 15). The largest relative error is 5.7% at Re = 15 and De = 1.
Re	[15, ∞]		[15, 4000]
Ca^a	[5, ∞]		[5, 500]
De	[0, ∞]^c		N/A^b
Mc	[0, 0.35]		[0.0064, 0.0068]
We	[13, ∞]		[100, 350]
	[0, 0.575]		[0.1, 0.25]

---

As discussed in Section 3, when the synthetically generated collapse time for n = 36 combinations of R_max ∈ [100, 400] μm and Λ_max ∈ [5, 9] are considered, pIMR identifies the Kelvin–Voigt model to be the optimal choice and recovers G and μ to within an accuracy of 5%. To examine the effect of the sample size n, we alternatively consider subsets with n = 9 and n = 3, as illustrated in Fig. 7(a). The corresponding results from pIMR are shown in Table 7. The n = 9 case results in Kelvin–Voigt parameters closely matching the n = 36 case, with the cost function decrement Δψ₀ decreasing from 2.66 to 0.20 when the constitutive model is advanced from Kelvin–Voigt to SLS. In contrast, the n = 3 case led to calibrated Kelvin–Voigt parameters with relative errors of 7% and 37%, respectively for G and μ. The minimized cost function ψ₀ decreases sharply from −5.06 to −31.18 when the relaxation time scale τ₁ is considered. This is due to the fact that exactly three experiments are considered to calibrate the three-parameter SLS model. Perhaps, a different subset of synthetic experiments with n = 3 would have resulted in a more accurate calibration of the viscoelastic model. However, such optimization of {R_max, Λ_max} is not feasible for real LIC experiments. As a general guideline, a large sample size of LIC experiments is beneficial for the performance of pIMR.

Fig. 7 Distribution of {R_max, Λ_max} in (a) synthetic experiments and (b) LIC experiments. For synthetic experiments, all data points are considered in the n = 36 case, the red square and bubble diamond points are considered in the n = 9 case, and only the blue diamond points are considered in the n = 3 case.

Table 7 Calibrated viscoelastic parameters, minimized cost function, and cost function decrement from synthetic experiments with varying levels of artificial errors. The input Kelvin–Voigt model parameters are {G = 10 kPa, μ = 0.10 Pa s}

Artificial error	n	Model	G (kPa)	μ (Pa s)	τ1 (μs)	ψ0	Δψ0
0%	36	NH	5.07	—	—	−2.79	—
		Newtonian	—	∼0	—	−1.94	—
		KV	10.11	0.095	—	−5.50	2.71
		SLS	10.38	0.116	0.506	−5.78	0.28
	9	NH	4.65	—	—	−2.76	—
		Newtonian	—	∼0	—	−1.99	—
		KV	10.11	0.095	—	−5.42	2.66
		SLS	10.10	0.104	3.34 × 10^−3	−5.62	0.20
	3	NH	5.46	—	—	−2.97	—
		Newtonian	—	∼0	—	−1.91	—
		KV	9.29	0.063	—	−5.06	2.09
		SLS	10.15	0.107	1.57	−31.18	26.12
+0.1%	36	NH	4.97	—	—	−2.79	—
		Newtonian	—	∼0	—	−1.95	—
		KV	10.00	0.094	—	−5.48	2.69
		SLS	10.27	0.116	0.515	−5.75	0.27
	9	NH	4.55	—	—	−2.76	—
		Newtonian	—	∼0	—	−2.01	—
		KV	10.00	0.095	—	−5.39	2.71
		SLS	10.34	0.118	0.603	−5.78	0.39
	3	NH	5.36	—	—	−2.96	—
		Newtonian	—	∼0	—	−1.92	—
		KV	9.18	0.063	—	−5.05	2.09
		SLS	10.05	0.107	1.583	−31.78	26.73
+1%	36	NH	4.07	—	—	−2.79	—
		Newtonian	—	∼0	—	−2.09	—
		KV	9.02	0.094	—	−5.27	2.48
		SLS	9.34	0.116	0.596	−5.48	0.21
	9	NH	3.65	—	—	−2.76	—
		Newtonian	—	∼0	—	−2.15	—
		KV	9.02	0.094	—	−5.21	2.45
		SLS	9.01	0.102	7.47 × 10^−4	−5.35	0.14
	3	NH	4.45	—	—	−2.97	—
		Newtonian	—	∼0	—	−2.05	—
		KV	8.18	0.062	—	−4.97	2.00
		SLS	9.12	0.109	1.686	−31.78	26.81

---

In our LIC experiment, the measurement of collapse time has a relative accuracy on the order of 0.1%. To examine the effect of such measurement uncertainty on pIMR, we repeat the above analysis with a relative error of 0.1% uniformly added to the collapse time of each experiment. Overall, the accuracy of the calibrated model parameters suffered minimally from the artificial error. In fact, for the n = 36 and n = 9 cases, the artificially increased collapse time led to a decreased G in the pIMR solution, matching the input value better than in the earlier, error-free case. When the artificial error is further increased to 1%, we observe that the calibrated G is decreased by approximately 11% compared to the error-free case, while the accuracy of μ shifted by less than 2%. This suggests that pIMR performs stably when processing collapse time data from our LIC experiments.

Acknowledgements

SHB and JBE acknowledge support from the U.S. Department of Defense, the Army Research Office under Grant No. W911NF-23-10324 (PMs Drs. Denise Ford and Robert Martin). MRJ acknowledges support from the U.S. Department of Defense under the DEPSCoR program Award No. FA9550-23-1-0485 (PM Dr Timothy Bentley). JBE, MRJ, and JY gratefully acknowledge support from the U.S. National Science Foundation (NSF) under Grant No. 2232426, 2232427, and 2232428, respectively.

References

1. J. B. Estrada, H. C. Cramer III, M. T. Scimone, S. Buyukozturk and C. Franck, Brain Multiphys., 2021, 2, 100034.
2. R. M. Wright, A. Post, B. Hoshizaki and K. T. Ramesh, J. Neurotrauma, 2013, 30, 102–118.
3. B. A. Abeid, M. L. Fabiilli, M. Aliabouzar and J. B. Estrada, J. Mech. Behav. Biomed. Mater., 2024, 160, 106776.
4. A. Vogel, S. Busch and U. Parlitz, J. Acoust. Soc. Am., 1996, 100, 148–165.
5. M. R. Bailey, R. O. Cleveland, T. Colonius, L. A. Crum, A. P. Evan, J. E. Lingeman, J. A. McAteer, O. A. Sapozhnikov and J. Williams, IEEE Symposium on Ultrasonics, 2003, vol. 2003, pp. 724–727.
6. Z. Xu, J. B. Fowlkes, E. D. Rothman, A. M. Levin and C. A. Cain, J. Acoust. Soc. Am., 2005, 117, 424–435.
7. J. E. Parsons, C. A. Cain, G. D. Abrams and J. B. Fowlkes, Ultrasound Med. Biol., 2006, 32, 115–129.
8. J. Vidal-Jove, X. Serres, E. Vlaisavljevich, J. Cannata, A. Duryea, R. Miller, X. Merino, M. Velat, Y. Kam and R. Bolduan, et al., Int. J. Hyperthermia, 2022, 39, 1115–1123.
9. M. Mendiratta-Lala, P. Wiggermann, M. Pech, X. Serres-Créixams, S. B. White, C. Davis, O. Ahmed, N. D. Parikh, M. Planert and M. Thormann, et al., Radiology, 2024, 312, e233051.
10. J. R. Sukovich, J. J. Macoskey, J. E. Lundt, T. I. Gerhardson, T. L. Hall and Z. Xu, IEEE Trans. Ultrason. Eng., 2020, 67, 1178–1191.
11. S. C. Haskell, E. Yeats, J. Shi, T. Hall, J. B. Fowlkes, Z. Xu and J. R. Sukovich, Ultrasound Med. Biol., 2025, 51, 909–920.
12. W. W. Chen, J. Dyn. Behav. Mater., 2016, 2, 2–14.
13. J. A. Zimberlin, N. Sanabria-DeLong, G. N. Tew and A. J. Crosby, Soft Matter, 2007, 3, 763–767.
14. S. Kundu and A. J. Crosby, Soft Matter, 2009, 5, 3963–3968.
15. M. P. Milner and S. B. Hutchens, Extreme Mech. Lett., 2019, 28, 69–75.
16. M. P. Milner and S. B. Hutchens, Mech. Mater., 2021, 154, 103741.
17. S. Raayai-Ardakani and T. Cohen, Extreme Mech. Lett., 2019, 31, 100536.
18. S. Chockalingam, C. Roth, T. Henzel and T. Cohen, J. Mech. Phys. Solids, 2020, 104172.
19. J. B. Estrada, C. Barajas, D. L. Henann, E. Johnsen and C. Franck, J. Mech. Phys. Solids, 2018, 112, 291–317.
20. J. Yang, H. C. Cramer III and C. Franck, Extreme Mech. Lett., 2020, 39, 100839.
21. J.-S. Spratt, M. Rodriguez, K. Schmidmayer, S. H. Bryngelson, J. Yang, C. Franck and T. Colonius, J. Mech. Phys. Solids, 2021, 152, 104455.
22. A. McGhee, J. Yang, E. Bremer, Z. Xu, H. Cramer III, J. Estrada, D. Henann and C. Franck, Exp. Mech., 2023, 63, 63–78.
23. J. B. Keller and M. Miksis, J. Acoust. Soc. Am., 1980, 68, 628–633.
24. A. Prosperetti and A. Lezzi, J. Fluid Mech., 1986, 168, 457–478.
25. R. Nigmatulin, N. Khabeev and F. Nagiev, Int. J. Heat Mass Transfer, 1981, 24, 1033–1044.
26. A. Prosperetti, L. A. Crum and K. W. Commander, J. Acoust. Soc. Am., 1988, 83, 502–514.
27. A. Prosperetti, J. Fluid Mech., 1991, 222, 587–616.
28. I. Akhatov, O. Lindau, A. Topolnikov, R. Mettin, N. Vakhitova and W. Lauterborn, Phys. Fluids, 2001, 13, 2805–2819.
29. C. Barajas and E. Johnsen, J. Acoust. Soc. Am., 2017, 141, 908–918.
30. X. Yang and C. C. Church, J. Acoust. Soc. Am., 2005, 118, 3595–3606.
31. C. Hua and E. Johnsen, Phys. Fluids, 2013, 25, 083101.
32. M. Warnez and E. Johnsen, Phys. Fluids, 2015, 27, 063103.
33. R. Gaudron, M. Warnez and E. Johnsen, J. Fluid Mech., 2015, 766, 54–75.
34. A. Sieber, D. Preso and M. Farhat, Phys. Fluids, 2023, 35, 027101.
35. J. L. Marsh, L. Zinnel and S. A. Bentil, J. Dyn. Behav. Mater., 2024, 1–17.
36. C.-D. Ohl, M. Arora, R. Ikink, N. De Jong, M. Versluis, M. Delius and D. Lohse, Biophys. J., 2006, 91, 4285–4295.
37. J. Yang, H. C. Cramer III, E. C. Bremer, S. Buyukozturk, Y. Yin and C. Franck, Extreme Mech. Lett., 2022, 51, 101572.
38. L. Mancia, E. Vlaisavljevich, N. Yousefi, M. Rodriguez, T. J. Ziemlewicz, F. T. Lee, D. Henann, C. Franck, Z. Xu and E. Johnsen, Phys. Med. Biol., 2019, 64, 225001.
39. J. B. Keller and I. I. Kolodner, J. Appl. Phys., 1956, 27, 1152–1161.
40. A. Tzoumaka, J. Yang, S. Buyukozturk, C. Franck and D. L. Henann, Soft Matter, 2023, 19, 3895–3909.
41. E. Bremer-Sai, J. Yang, A. McGhee and C. Franck, Exp. Mech., 2024, 64, 587–592.
42. J. W. Strutt, London, Edinburgh Dublin Philos. Mag. J. Sci., 1917, 34, 94–98.
43. J. Yang, A. McGhee, G. Radtke, M. Rodriguez and C. Franck, Phys. Fluids, 2024, 36, 017136.
44. C. E. Brennen, Cavitation and bubble dynamics, Cambridge University Press, 2014.
45. C. Zener, Elasticity and anelasticity of metals, University of Chicago Press, 1965.
46. J. R. Sukovich, S. C. Haskell, Z. Xu and T. L. Hall, J. Acoust. Soc. Am., 2020, 147, 1339–1343.
47. J. R. Tse and A. J. Engler, Curr. Protoc. Cell Biol., 2010, 47, 10–16.
48. S. Buyukozturk, J.-S. Spratt, D. Henann, T. Colonius and C. Franck, Exp. Mech., 2022, 62, 1037–1050.
49. J. A. Nelder and R. Mead, Comput. J., 1965, 7, 308–313.
50. J. C. Lagarias, J. A. Reeds, M. H. Wright and P. E. Wright, SIAM J. Optim., 1998, 9, 112–147.
51. Z. Wang, X. Huan and K. Garikipati, arXiv, 2018, preprint, arXiv:1812.11285 DOI:10.48550/arXiv.1812.11285.
52. Z. Wang, J. B. Estrada, E. M. Arruda and K. Garikipati, J. Mech. Phys. Solids, 2021, 153, 104474.
53. D. P. Nikolov, S. Srivastava, B. A. Abeid, U. M. Scheven, E. M. Arruda, K. Garikipati and J. B. Estrada, Philos. Trans. R. Soc., A, 2022, 380, 20210324.
54. R. Tibshirani, J. R. Stat. Soc. Ser. B: Stat. Methodol., 1996, 58, 267–288.
55. M. D. McKay, R. J. Beckman and W. J. Conover, Technometrics, 2000, 42, 55–61.
56. D. G. Strange and M. L. Oyen, J. Mech. Behav. Biomed. Mater., 2012, 11, 16–26.
57. W. Gu, H. Yao, C. Huang and H. Cheung, J. Biomech., 2003, 36, 593–598.
58. V. Normand, D. L. Lootens, E. Amici, K. P. Plucknett and P. Aymard, Biomacromolecules, 2000, 1, 730–738.
59. K. Takamura, H. Fischer and N. R. Morrow, J. Pet. Sci. Eng., 2012, 98, 50–60.
60. T. Chu, J. B. Estrada and S. H. Bryngelson, Comput. Mech., 2025, 1–17.
61. T. Chu, J. B. Estrada and S. H. Bryngelson, arXiv, 2025, preprint, arXiv:2505.13283 DOI:10.48550/arXiv.2505.13283.
62. S. Bentil, K. Ramesh and T. Nguyen, Exp. Mech., 2016, 56, 759–769.
63. R. S. Salzar, D. Treichler, A. Wardlaw, G. Weiss and J. Goeller, J. Neurotrauma, 2017, 34, 1589–1602.
64. H. Saraf, K. Ramesh, A. Lennon, A. Merkle and J. Roberts, Exp. Mech., 2007, 47, 439–449.
65. J. R. Sukovich, Z. Xu, Y. Kim, H. Cao, T.-S. Nguyen, A. S. Pandey, T. L. Hall and C. A. Cain, IEEE Trans. Ultrason. Eng., 2016, 63, 671–682.
66. J. J. Macoskey, T. L. Hall, J. R. Sukovich, S. W. Choi, K. Ives, E. Johnsen, C. A. Cain and Z. Xu, IEEE Trans. Ultrason. Eng., 2018, 65, 2073–2085.
67. E. Traldi, M. Boselli, E. Simoncelli, A. Stancampiano, M. Gherardi, V. Colombo and G. S. Settles, EPJ Tech. Instrum., 2018, 5, 1–23.
68. V. Agrež, T. Požar and R. Petkovšek, Opt. Lett., 2020, 45, 1547–1550.

---

Footnote

† These authors contributed equally to this work.

---