Shell viscosity estimation of lipid-coated microbubbles

Understanding the shell rheology of ultrasound contrast agent microbubbles is vital for anticipating their bioeffects in clinical practice. Past studies using sophisticated acoustic and optical techniques have made enormous progress in this direction, enabling the development of shell models that adequately reproduce the nonlinear behaviour of the coated microbubble under acoustic excitation. However, there have also been puzzling discrepancies and missing physical explanations for the dependency of shell viscosity on the equilibrium bubble radius, which demands further experimental investigations. In this study, we aim to unravel the cause of such behaviour by performing a refined characterisation of the shell viscosity. We use ultra-high-speed microscopy imaging, optical trapping and wide-field fluorescence to accurately record the individual microbubble response upon ultrasound driving across a range of bubble sizes. An advanced model of bubble dynamics is validated and employed to infer the shell viscosity of single bubbles from their radial time evolution. The resulting values reveal a prominent variability of the shell viscosity of about an order of magnitude and no dependency on the bubble size, which is contrary to previous studies. We find that the method called bubble spectroscopy, which has been used extensively in the past to determine the shell viscosity, is highly sensitive to methodology inaccuracies, and we demonstrate through analytical arguments that the previously reported unphysical trends are an artifact of these biases. We also show the importance of correct bubble sizing, as errors in this aspect can also lead to unphysical trends in shell viscosity, when estimated through a nonlinear fitting from the time response of the bubble.


I. INTRODUCTION
During the past decades, the development of stabilised gas microbubbles smaller than a red blood cell (< 10 µm), and thus able to navigate the thinnest capillaries, has opened unique opportunities to utilise the effects of ultrasonic cavitation in medicine. Upon high-amplitude ultrasound driving, the microbubbles are set into a periodic succession of growth and collapse phases. The microbubbles thus pump energy out of the acoustic wave and release it locally in the form of sound (even at frequencies other than the driving frequency) and mechanical action. This extraordinary energy-focusing ability makes a simple bubble an extremely capable device for a large number of biomedical applications. Coated microbubbles are currently used as ultrasound contrast agents to improve the contrast in perfusion imaging of the myocardium [1,2], liver [3], brain [4], kidney [5] and other organs, providing unique information on blood volume and velocity owing to their large echogenicity. The imaging potential of microbubbles can be augmented by equipping the coating with specific ligands. Site-targeted contrast agents can be used to non-invasively image molecular events in vivo such as inflammation [6,7], angiogenesis [8,9] and early tumour formation [10,11]. Moreover, microbubbles are presently investigated for therapeutic applications. The high-energy mechanical effects of bubble cavitation (e.g. acoustic microstreaming [12] and microjetting [13]) can be harnessed to promote clot lysis in acute ischemic stroke [14,15], deliver plasmid DNA and drugs in the immediate perivascular region for myocardial infarction [16], atherosclerosis [17] and tumours treatment [18,19] or to locally permeabilise the blood-brain barrier to enhance the delivery of therapeutic agents for brain cancer [20][21][22], Alzheimer's disease [23] and amyotrophic lateral sclerosis [24] treatment.
A prolonged intravascular stability against dissolution and coalescence is achieved by encapsulating the microbubbles with a protein, biopolymer or lipid shell and by using a high-molecular-mass gas content such as perfluorocarbons and sulphur hexafluoride, which are less prone to outward loss compared to air [25]. To date, most clinically approved ultrasound contrast agents are coated with a phospholipid shell which is more flexible than other types of coatings, yielding less stable but more echogenic microbubbles. Lipids naturally adsorb on gas-liquid interfaces and self-assemble into a monolayer which lowers the gas-liquid surface tension and thus the large Laplace pressure, arresting the gas efflux. Moreover, it imparts rheological properties to the surface [26]. Monolayers of DPPC (1,2-distearoyl-sn-glycero-3phosphocholine), a frequent main lipid constituent of the bubble shell, show at low harmonic driving frequencies a dominantly viscous behaviour similar to threedimensional liquids. On the other hand, at higher frequencies their microstructure dominates the interfacial rheology and induces an additional elastic contribution [27]. This elasticity is of intrinsic rheological origin and arXiv:2302.06339v2 [physics.flu-dyn] 11 Aug 2023 is added to the Gibbs (apparent) monolayer elasticity, which arises when interface deformations yield variations of surfactant concentration [28]. The properties of the shell have a fundamental impact on the dynamical behaviour of ultrasound contrast agent microbubbles, and consequently, on their biophysical effects. Compared to uncoated bubbles, coated bubbles show a higher resonance frequency due to a stiffer interface, a lower amplitude of response owing to viscous dissipation mechanisms between lipid molecules and a plethora of nonlinear effects [29].
Unfortunately, characterising the shell of a microbubble under clinically-relevant working conditions poses significant challenges owing to the micrometric size and the extremely high driving frequencies (∼ 10 6 Hz) to which these bubbles are subjected. Direct experimental rheological techniques such as Langmuir trough and capillary pressure tensiometry methods are limited to significantly lower driving frequencies (∼ 1 Hz and ∼ 10 2 Hz, respectively) [30] and at least one order of magnitude smaller film curvatures [31]. Therefore, as current direct approaches are not adequate for characterising the shell properties, only indirect methods can be used. The seminal works of de Jong et al. [32], Frinking & de Jong [33], Hoff et al. [34] and Gorce et al. [35] provided the first estimates of the elastic and viscous contribution of the shell by performing acoustic attenuation measurements of polydisperse suspensions of coated microbubbles. This approach was later refined by Parrales et al. [36] and Segers et al. [37] by employing monodisperse bubble suspensions produced by microfluidic flow-focusing techniques. Moreover, Helfield et al. [38] estimated the shell properties from scattered pressure measurements of single bubbles. Morgan et al. [39] provided the first optical characterisation by fitting a shell model to the radial evolution, captured on a streak camera image, of single ultrasound-driven microbubbles. The same system was later also used by Doinikov et al. [40]. Only with the development of the groundbreaking Brandaris camera [41] was the curtain on the wide-field time-resolved dynamics of microbubble finally lifted. van der Meer et al. [42] and Overvelde et al. [43] were the first to leverage this new potential to characterise the shell of single DSPCbased microbubbles, followed by Luan et al. [44] to investigate the effect of liposomes loaded on the shell and van Rooij et al. [45] to evaluate the influence of the carbon chain length of the phospholipids on the shell properties. A side laser scattering system was used by Tu et al. [46,47] to characterise the shell properties of commercial phospholipid-coated microbubbles by measuring their response to an ultrasound pulse. The sensitivity of this technique has been further improved by employing a forward laser-scattering system that enabled Lum et al. [48] to measure the decaying free oscillations and Lum et al. [49], Supponen et al. [50] and Daeichin et al. [51] to retrieve the frequency response of single bubbles set into small-amplitude oscillations via photoacoustic forcing. Such small oscillations (few tens of nanometers) simplify the characterisation, as the absence of nonlinear effects allows the use of linear models.
In contrast to the dilatational modulus, for which there is broad empirical evidence regarding its independence from the bubble size, the reported dilatational viscosity values have been shown to increase with the bubble equilibrium radius and with a growing rate that strongly varies from study to study [36-40, 42, 44-47, 50, 51]. This observed dependency is widely regarded as unphysical, as it is illogical for a material property to be dependent on the quantity of the material. van der Meer et al. [42], Doinikov et al. [40] and Tu et al. [47] also showed the shell viscosity to decrease with the (maximum) strain rate of the bubble, suggesting a strainthinning behaviour to possibly explain its questionable dependency on the bubble size. However, this conclusion contradicts the observations from Helfield et al. [38] and van Rooij et al. [45] who found no significant relationship between shell viscosity and strain rate. Moreover, the results reported by Tu et al. [47] are particularly puzzling: as the parameters constituting the maximum strain rate (Ṙ/R) max ≈ 2πf ∆R max /R 0 -bubble resting radius R 0 , maximum radial expansion (through pressure) ∆R max and driving frequency f -vary, the shell viscosity gives rise to different families of curves instead of collapsing on the same universal curve. In our view, this means that the dependence of shell viscosity on the initial radius cannot be explained by a dependence on the strain rate. Furthermore, the assumption that the shell viscosity is dependent on the strain rate throughout the entire radial excursion of the bubble pursuant to the Cross model, as proposed by Doinikov et al. [40] following the result from van der Meer et al., is questionable. In fact, during the expansion phase of the bubble, the phospholipids separate from each other, leading to a loss of the microstructure responsible for the non-Newtonian characteristics of the shell beyond a certain level of interfacial strain. Conversely, during the compression phase of the bubble, the in-plane strain of the microstructure may have a minor role compared to the out-of-plane buckling. Most studies reporting a dependency of shell viscosity on the bubble size have estimated the shell properties by analysing the bubble response in the frequency domain, adopting the "bubble spectroscopy" technique introduced by van der Meer et al. [42] or a variation of it. The remaining works instead inferred the shell properties by fitting the experimental bubble radius-time curves.
In this work, we aim to clarify the apparent dependency of the shell viscosity on the microbubble size and the unexplained discrepancies in the values reported in different studies for similar bubbles. We use ultra-highspeed microscopic imaging and optical trapping to directly observe the response of single ultrasonically driven coated microbubbles in an unbound fluid environment across a range of bubble radii. We characterise the shell dilatational viscosity examining the bubble response in the time domain, instead of adopting the common bubble spectroscopy approach. Contrary to other studies which fitted the bubble time response, we first confronted our experimental data against the theoretical predictions to verify that the chosen bubble model is capable of describing the salient features of the bubble response, in particular the bubble resonance size. The shell dilatational viscosity is then inferred by fitting the experimental response with the theoretical one. Finally, the results are compared with previous studies and an explanation for the origin of the unphysical dependence of shell viscosity on the microbubble size is provided.

II. THEORY
The radial motion of a spherical bubble of radius R(t) is described by a modified form of the Rayleigh-Plesset equation for mildly compressible Newtonian media [52][53][54][55], which includes a pressure term Σ R,Ṙ to account for the generalised interfacial stresses σ s R,Ṙ and reads where over-dots denote time differentiation, ρ l the liquid density, c l the speed of sound in the medium, p g the gas pressure inside the bubble, p ∞ the undisturbed ambient pressure, p d (t) the acoustic driving pressure and µ l the dynamic viscosity of the medium. The pressure contribution of the interface Σ R,Ṙ is given by the radial component of the surface divergence of the interfacial stress [26,29] where ∇ s def = I s · ∇ is the surface gradient operator, I s def = I − nn is the surface unit tensor and n is the unit vector normal to the interface. As opposed to earlier models that treated the shell as a three-dimensional solid [34,39,56], Chatterjee and Sarkar [57,58] were the first to suggest modelling the bubble coating as a twodimensional continuum. This approach is considered the most appropriate for molecular membranes (i.e. phospholipid [59] and protein [60] shells) for which the assumption of continuum in the normal direction is questionable [61]. They proposed to model the microbubble shell as a viscoelastic solid interface, expressing the surface stress tensor, owing to the spherical symmetry of the problem, as The first term represents the isotropic part of the Evans-Skalak model for elastic interfaces [62] without differentiating the thermodynamic contribution from the rheological one [63,64] and the second the isotropic part of the Boussinesq-Scriven model for viscous interfaces [26] which is the surface equivalent of the Newtonian model for compressible bulk fluids. σ 0 is the interfacial surface tension at the equilibrium bubble radius R 0 , E s is the interfacial dilatational modulus, κ s is the interfacial dilatational viscosity and J = R 2 /R 2 0 is the relative area deformation. Upon substitution of Eq. (3) into Eq. (2) and by making use of the following differential identities, valid for a spherically symmetric problem, ∇ s · n = 2/R and ∇ s · I s = − (2/R) n, the pressure contribution given by the interface reads To account for the strong nonlinear effects that coated microbubbles exhibit, such as subharmonic oscillations [65][66][67][68][69], compression-only behaviour [70], thresholding behaviour [71] and resonance frequency shifting [43], which a quasi-linear shell model as Eq. (4) cannot represent, Marmottant et al. [70] proposed a phenomenological extension inspired by the behaviour of liquid-gas interfaces covered by insoluble surfactants [27], yielding which accounts for: (i) the buckling of the coating when the bubble is compressed below its shell compression limit radius R buckling resulting in a tensionless interface, (ii) the rupturing of the coating when the bubble is expanded beyond R rupture where it experiences a surface tension equal to that of a pure gas-water interface σ water = 72.8 mN m −1 , and (iii) the elastic regime in between. The Marmottant model has proven to be a compelling engineering representation of the microbubble shell behaviour, capable of explaining a broad spectrum of nonlinearities [43,69], and to be sophisticated enough to accurately predict the behaviour of bubbles in the present study, without the need for additional free parameters. However, this model is not entirely free from physical inconsistencies, an important one being the assumption of a constant shell viscosity for varying the free molecular area. Nevertheless, this can be justified as being an average viscosity over the radial excursion of the bubble.
In order to close the problem, the gas pressure contained in the bubble must be determined, in principle, by solving the full Navier-Stokes equation inside the bubble. This can be done analytically only for small-amplitude oscillations for which the governing equations can be linearised [72,73] or numerically at a high computational cost. For these reasons, extensive use is made of the crude approximation of a polytropic process: with p g,0 the bubble internal pressure at rest and n the polytropic index. In spite of its appealing simplicity, this approximation suffers from serious limitations [74]. First, its range of validity is restricted to low-amplitude oscillations and to global Péclet numbers Pe = R 2 0 ω/D g (where ω = 2πf is the driving angular frequency and D g is the gas thermal diffusivity) much smaller or much bigger than one, and therefore to the particular cases where the gas follows an isothermal process (n = 1) or an adiabatic one (n = γ, where γ is the gas specific heat ratio), respectively, for the whole radial motion of the bubble. And second, its use results in no energy loss associated with the change of temperature of the gas. These restrictions are unfortunate in the context of ultrasound contrast agent microbubbles because under clinical conditions they often undergo large nonlinear oscillations, characterised by a thermal motion defined by a global Pe ∼ 1 with both isothermal and adiabatic behaviours occurring at different phases of oscillation, and subjected to a thermal loss of the same order of magnitude of the viscous damping. To address the latter point, it is customary in the literature to artificially increase the liquid viscosity doubling its value [75]. These considerations motivate us to take a step beyond the polytropic approximation. In an attempt to simplify the problem, Nigmatulin et al. [76] and Prosperetti et al. [74] independently found the spatial variation of the pressure inside the bubble to be negligible in most cases (∆p g /p g ∼ O(RṘ/λc g ,Ṙ 2 /c 2 g ), where λ and c g are representative values of the sound wavelength and speed in the bubble), allowing, for a perfect gas, to exactly express the gas radial velocity as: and from this result recover an exact relationship for the gas pressure: where r is the radial coordinate, T g is the gas temperature, and K g is the thermal conductivity of the gas. These two relationships allow to reduce the number of governing partial differential equations (PDEs) from three to one, the one necessary to find the temperature gradient at the bubble surface, i.e. ∂ r T g | R , which is the term that governs the heat exchange and thus closes the problem. One can solve either the energy equation or the continuity equation coupled with the equation of state for perfect gases. The first alternative was chosen by Kamath & Prosperetti [77], who solved the energy equation by using a Galerkin-Chebyshev spectral method. Recently, the second alternative was pursued by Zhou [78], who found the temperature profile inside the bubble to be divided into three regions: (1) an internal layer characterised by uniform temperature, (2) a buffer layer, and (3) a layer adjacent to the bubble surface characterised by a linear temperature distribution, which in turn suggests dividing the computational domain into as many cells. In this way, the continuity equatioṅ where ρ g denotes the gas density, can be converted into two ordinary differential equations (ODEs) where m g,i is the gas mass in cell i and f j is the mass flux across the cell interface j. ρ uw g,j is the density of the neighboring cell of interface j in the upwind direction, u rel g,j is the convective velocity and S j is the interface j surface area. The velocity u g can be obtained from Eq. (7), while the temperature T g computed through the equation of state for ideal gas where R is the the specific gas constant. The change in bubble surface temperature is only a small fraction of that of the gas and can, therefore, be neglected, i.e. T g | R ≈ T l [74]. At this point, the motion of the bubble can be effortlessly solved by integrating in time the system of ODEs consisting of Eq. (1), (8) and (10). The reader can refer to the article by Zhou [78] for full details of the method. The simplicity and accuracy of the model make it suitable for this work.

A. Microbubble synthesis
Lipid-coated microbubbles with a C 4 F 10 (perfluorobutane, Fluoromed) gas core are prepared in-house by probe sonication. The lipid coating consists of 90 mol% DPPC (1,2-distearoyl-sn-glycero-3-phosphocholine, NOF EUROPE) and 10 mol% of DPPE-PEG5K (1,2-dipalmitoyl-sn-glycero-3-phosphoethanolamine-N-[methoxy(polyethylene glycol)-5000]), Avanti Polar Lipids). The lipids are first dissolved in chloroform. The solution is placed under vacuum at 35°C overnight to remove the solvent and produce a dry lipid film. The dry lipid film is rehydrated with 1× PBS (phosphate buffered saline, Boston BioProducts) to yield a total lipid concentration of 2 mg ml −1 . The solution is then sonicated (SFX550, Branson; 20 kHz, 550 W) at low power (30%) for five minutes to convert the multilamellar vesicles to unilamellar liposomes. Microbubbles are formed by probe-sonicating the surface of the lipid solution at full power for ten seconds while simultaneously flowing C 4 F 10 gas over it. The microbubble suspension is then cooled down to room temperature and washed using centrifugation. Differential centrifugation is used to size-isolate the microbubbles into a specific size distribution (1 − 4 µm-radius), based on the protocol developed by Feshitan et al. [79]. Fluorescent-labeled microbubbles are fabricated by adding the lipid dye DiI (1,1'-dioctadecyl-3,3,3'3'-tetramethylindocarbocyanine perchlorate, Sigma-Aldrich) to the lipid solution at a concentration of 5 µg ml −1 before probe sonication.

B. Experimental setup
To study the dynamics of microbubbles, we developed a multimodal microscope ("Mscope", Fig. 1) which allows for wide-field ultra-high-speed imaging, holographic optical trapping and wide-field fluorescence.
The main body of the system [ Fig. 1(M)] is a custombuilt upright microscope realised using modular optomechanics components (Thorlabs, cage system), installed on an optical table (T1220C, Thorlabs) with active isolators (PTS603, Thorlabs) to minimise environmental vibrations. The microscope is equipped with a waterdipping objective lens of focal length F = 2 mm (CFI Plan 100XC W, Nikon) and a F = 400 mm tube lens (TL400-A, Thorlabs) for a total magnification of 200×. An ultra-high-speed camera (HPV-X2, Shimadzu) allows for recordings at 10 million frames per second over 25.6 µs of continuous visualisation of a 64 µm × 40 µm field of view with a 160-nm pixel resolution. Backlight illumination is provided by a continuous halogen illuminator (OSL2, Thorlabs) for live imaging and two Xenon flash lamps used sequentially (MVS-7010, EG&G) for video recording, all combined into a single optical fiber output and focused on the sample through a custom-built condenser (L1 and L2, F = 25 mm, AC127-025-A, Thorlabs). A water bath (T l ≈ 22°C) filled with deionised water is accommodated at the base of the microscope. A custom-designed 3D-printed test chamber with optically and acoustically transparent windows on top and bottom (low-density polyethylene, acoustic impedance z = 1.79 × 10 6 Pa s m −1 ) is suspended in the water bath and its position controlled by a three-axis motorised microtranslation stage (PT3/M-Z8, Thorlabs). A broadband focused ultrasound transducer (PA1612, Precision Acoustics; 1.5 MHz of center frequency, 85 mm of focal length, ≈ 4 mm of beam width at −6 dB) is positioned in the water bath at an angle of 30 • with respect to the horizontal plane to minimise acoustic reflections and manoeuvred using a manual three-axis microtranslation stage (#12694, #66-511, Edmund Optics). A function generator (LW 420B, Teledyne LeCroy) is used to produce the driving pulse, which is then amplified by a radiofrequency power amplifier (1020L, E&I). A calibrated needle hydrophone (0.2 mm, NH0200, Precision Acoustics) is employed to measure the driving acoustic pressure inside the test chamber and to align the acoustic focal point with the optical field.
The holographic optical tweezers and fluorescent microscopy are realised using a 532-nm optically pumped semiconductor continuous-wave laser (Verdi G10, Coherent). The laser beam enters an acousto-optic tunable fil- The transmitted beam (0 th order) is ceased with a beam dump while the diffracted beam (1 st order), whose intensity is dependent upon the radiofrequency drive power provided by a radiofrequency driver (MOD.8C.10.b.VIS, AA Opto Electronic), is selectively sent to the optical tweezers or the fluorescence setup by means of a movable mirror.
Three-dimensional optical trapping of single microbubbles has been demonstrated in the pioneering works of Prentice et al. [80] and Garbin et al. [81] by focus-ing an optical vortex laser beam obtained by converting the laser transverse mode from Gaussian to Laguerre-Gaussian by means of a diffractive optical element. Compared to the infrared laser employed in those studies, using a green laser reduces the local heating of water by three orders of magnitude [82,83]. In our implementation [ Fig. 1(HOT)], the optical path of the optical tweezers setup is kept as short as possible to achieve the highest pointing stability. The diffractive optical element (DOE) is provided by a reflective spatial light modulator (SLM) (PLUTO-2.1 VIS-096, HOLOEYE Photonics). The SLM is placed after a 7.5× galileian telescope, made of a pair of two lenses (L3, F = −20 mm, ACN127-020-A, Thorlabs; L4, F = 150 mm, AC254-150-A, Thorlabs), and an iris (SM1D12, Thorlabs) which expand and spatially filter the laser beam to completely fill the active area of the SLM with an approximately uniform-intensity distribution which helps to generate a strong optical trap. The linear polarisation of the laser is rotated by using a zero-order half-wave plate (WPHSM05-532, Thorlabs) to match the alignment of the liquid crystalline molecules of the SLM microdisplay. The laser beam mode is converted into a Laguerre-Gaussian one by implementing an optical vortex phase mask on the SLM. In order to minimise the diffraction effects, a 4F system is used to relay the SLM image plane to the objective back focal plane (L5, F = 400 mm, AC254-400-A, Thorlabs; L6, F = 200 mm, AC254-200-A, Thorlabs). The magnification of the 4F system is 0.5× in order to slightly overfill the objective lens back aperture and thus obtain the best optical trapping performance. The 0 th -order component of the beam, unavoidably present due to the pixelated structure of the SLM, is suppressed by performing spatial filtering after having separated the 1 st -order-component image plane from the 0 th one by superimposing a phase mask for a lens onto the SLM. Upon reflection on a reflective-band dichroic beamsplitter (ZT532dcrb, Chroma), the beam is focused on the sample by the objective lens. The optical trap can be freely positioned within the sample image plane by superimposing a phase mask for a prism onto the SLM.
With the use of movable mirrors the optical tweezers setup can be bypassed and make use of the laser for performing epifluorescence microscopy [ Fig. 1(F)]. The size of the field of illumination and consequently the laserexcitation power density can be adjusted varying the spacing of a lens relay system (L7 and L8, F = 100 mm, AC254-100-A, Thorlabs). Undesirable light wavelengths are prevented from illuminating the specimen by means of a narrow passband excitation filter (ZET532/10x, Chroma). The dichroic beamsplitter reflects the excitation light into the objective lens. The laser line is removed from the specimen image with an emission filter (ET590/50m, Chroma).

A. Microbubble radius measurement
Measuring accurately the radius of a microbubble is essential to equitably compare its experimental behaviour with theoretical prediction. However, making this measurement with bright-field microscopy is not trivial as the size of a microbubble is comparable to the wavelength of light. The light waves passing around the microbubble generate a Fresnel diffraction pattern that appears as concentric bright and dark circles around the bubble and thus hinders the accurate estimation of its dimension. Fortunately, fluorescence microscopy does not suffer from the interference between the light rays and the specimen because the specimen itself is selfemitting. By adding a fluorophore to the shell composition, it is possible to leverage fluorescence microscopy to image the microbubble. However, the brightness intensity of fluorophores is not sufficient for ultra-high speed recordings. Fluorescence microscopy is therefore only employed to evaluate the systematic error present in the bubble radius measurements obtained with bright-field microscopy. The error is assessed by comparing the bubble radius as measured with fluorescence and bright-field microscopy. Figure 2 compares the two, showing that bright-field microscopy overestimates the bubble radius by R err = (140 ± 40) nm across all bubble sizes. It is important to mention that this error value is characteristic to this specific optical setup and follows from the particular combination of the light source and objective lens numerical apertures used. It should not be surprising that a discrepancy smaller than the pixel resolution can be recovered, owing to the sub-pixel accuracy of the edge detection algorithm employed for detecting the bubble contour.

B. Microbubble response upon ultrasonic driving
The lipid-coated microbubbles are injected at a very low dilution ratio in the test chamber in order to avoid inter-bubble radiative forces (secondary Bjerknes forces) upon acoustic excitation. Because of buoyancy, the bubbles float up and adhere to the top window of the test chamber. Individual bubbles are then optically trapped and moved away from the top window by > 50 µm using optical tweezers in order to avoid wall interferences. The order of the optical vortex and the laser intensity are adjusted manually for each bubble size. The bubble under investigation is acoustically driven by a 20-cycle sinusoidal pulse at 1.5 MHz and with 40 kPa of peak negative pressure produced by the ultrasound transducer. A video recording at 10 million frames per second of the bubble response is performed using the ultra-high-speed videomicroscopy apparatus. The optical trap is deactivated during the recording by using the microsecond-fast electronic shutter. chain of events occurring during a single experiment, including the measured pressure driving pulse (black), the measured laser intensity (i.e., the optical trap strength) (green) and the recording time window (orange).
Radius-time curves are extracted from the recordings using a feature-extraction image-processing algorithm. The normalised maximum bubble radial expansion R max /R 0 − 1 is taken as the exemplary feature of the bubble response and as the target for the inverse problem of inferring the shell viscosity from the bubble response. This choice is more robust than considering the average maximum expansion over several oscillation cycles because it prevents the bubble dissolution, which may occur in the final cycles to a varying extent, or shape modes, which develop for large bubbles again in the final cycles, from influencing the bubble response evaluation and the consequent shell viscosity estimation. In fact, bubble shrinkage leads to a reduction in radial excursion while shape modes cause the loss of the bubble's spherical symmetry and thus the inability to correctly relate the bubble's cross sectional area with its volume and hence to an equivalent radius. Nevertheless, these issues always occur after the maximum bubble expansion is reached, which we can easily verify visually from individual videos. By selecting the maximum bubble expansion as the target for the inverse problem of inferring the shell viscosity from the bubble response, we effectively prevent these concerns from compromising the estimation of shell viscosity. Figure 3(b) depicts with orange dots the distribution of the normalised maximum bubble radial expansion across a range of bubble equilibrium radii R 0 . The bubble radius is adjusted by considering the systematic sizing error coming from the Fresnel interference and deemed equal to R err = 140 nm.
As expected, the curve exhibits a pronounced peak at a specific bubble radius that represents the resonance bubble radius. However, it also unveils a large variability in the response for microbubbles of similar size. A representative example of this arbitrariness is detailed in Fig.  3(c), which shows the large difference in the radial time evolution of two bubbles of nearly identical size (R 0 = 2.5 µm), reporting their radius-time curves in solid orange lines and selected snapshots of their dynamics, extracted from Supplementary Movies 1 and 2, in the orange boxes. This irregular bubble behaviour can be explained with a dispersion in the shell parameters: dilatational modulus but most prominently dilatational viscosity. The orangeshaded area in Fig. 3(b), computed using the bubble dynamics model reported in §II, shows how a range of dilatational modulus E s = 0.3 − 2 N m −1 and dilatational viscosity κ s = 1.5 × 10 −9 − 12 × 10 −9 kg s −1 allows us to fully describe the variability in the experimental results. This figure also suggests a useful qualitative classification of the bubble response into two main regimes, i.e., shelldominated and inertia-dominated, based on the bubble size.
In the shell-dominated regime, the shell properties have a large influence on the bubble dynamics. Therefore, minor changes in the shell rheology may lead to profoundly different bubble response outcomes. This leverage effect is amplified in the resonance sub-regime where a bubble experiences the largest radial excursion. It should be noted that high radial excursion values can possibly result in the elastic shell regime to play a marginal role in the overall bubble response compared to the prevailing buckling and rupture regimes. This is also the case for the combination of shell composition and negative peak pressure employed in this study. In fact, as shown in Fig. 3(b), the bubble response in the resonance regime is not influenced by the choice of the dilatational modulus. The grey-shaded area, representing the dispersion of the bubble response for the same range in shell viscosity reported before but now for a specific dilatational modulus value E s = 0.6 N m −1 , matches finely the orange-shaded one. Therefore, for clinical applications, where the bubble suspension is subjected to relatively strong ultrasound forcing and the resonance is leveraged to get the maximum result at the minimum mechanical index, characterising the shell viscosity is certainly more important than doing so for the dilatational modulus. Nonetheless, for sub-resonant bubbles, the response is less sensitive to shell viscosity variations and therefore these alone -even with a larger range of values than the one used for Fig. 3(b) -cannot fully explain the scatter found in the experimental results, which must therefore also be attributed to the dilatational modulus.
On the other hand, in the inertia-dominated regime, the shell plays a marginal role and the radial excursion of the bubble is mainly dependent on its equilibrium radius, and therefore on its inertia, as evident from the different curves representing different values of shell properties col- lapsing on a single one. As a consequence, the elevation of this curve can only be controlled by the driving pressure. Therefore, finding a good match between experiments and theory in this regime means that the measured and theoretical pressure are coinciding, propping our acoustic pressure measurement. Also note how the variability of the bubble response drops in this regime which correlates with the diminished influence of shell properties on the response itself. This observation further reinforces the idea that the observed variability in the bubble response stems from the scattering of shell properties from one bubble to another.
The measured bubble radial time evolutions agree well with the theoretical predictions as shown by the representative radius-time curves in Fig. 3(c). The bubble with the largest excursion is best modelled with a shell viscosity κ s = 1.8 × 10 −9 kg s −1 . The discrepancy in the final cycles and the new equilibrium radius are caused by the net gas loss of the bubble, which is likely facilitated by the partial shedding of the shell. These two phenomena are not accounted for in the theoretical model. Consequently, as the ultrasound burst progresses, the bubble experiences changes in its resonance frequency and shell properties compared to the initial conditions, leading to deviations in its response from the theoretical predictions. Conversely, the bubble with the smallest radial excursion is simulated with a much higher shell viscosity κ s = 10 × 10 −9 kg s −1 . In this case, the bubble does not undergo shrinkage, resulting in excellent agreement between the theoretical and experimental curves, even in the tail region. This result further strengthens our hypothesis regarding shell viscosity being the prime actor in explaining the variability in bubble response. We emphasise that bubble dissolution only occurs during intense radial oscillations. Otherwise, the bubbles remain stable for extended periods, tens of minutes, within the aqueous solution.
The theoretical bubble response is calculated by setting the initial surface tension σ 0 to zero because no finite value can encompass all data points, whereas a null value successfully does so, as illustrated in Fig.  3(b). It is worth noting that a tensionless or nearly tensionless (null Laplace pressure) bubble is also consistent with its long-term stability against dissolution in a saturated medium [25]. Although the microbubbles are originally fabricated with a perfluorobutane gas core, it is likely that this has been replaced by air due to the extended period of time (few tens of minutes) that the bubbles reside in an air-saturated aqueous medium. A study by Kwan and Borden [84], supports this notion, revealing that the heavy-molecular-weight gas initially present in the gas core of microbubbles is substituted by the air dissolved in the surrounding aqueous medium within a few minutes. Accordingly, we set γ = 1.

C. Shell dilatational viscosity estimation
The shell dilatational viscosity of each tested microbubble is inferred from the observation of its radial time evolution by solving an inverse problem. A grid search optimisation algorithm is used to find the theoretical bubble radial time evolution that results in the minimum deviation between empirical and theoretical maximum bubble expansions. The free parameters are the dilatational shell viscosity κ s and modulus E s . The other parameters are held fixed at the values presented in the previous section. We would like to emphasise that due to the typical shape of pressure signals generated by ultrasound transducers, inferring shell viscosity from the maximum expansion of the bubble is a more robust approach than inferring it from the decaying rate of the bubble response tail. The response tail, in fact, does not truly represent the free dissipative relaxation of the bubble but rather the forced response to the tail of the pressure driving signal [see Fig. 3(a)]. On the contrary, the (maximum) radial expansion of a bubble correlates robustly with its shell viscosity as shown in Fig. 3(b). Figure 4 reports with the orange dots the values of shell viscosity resulting from solving the inverse problem, which show a large variability of one order of magnitude and an average value of κ s = 5.4 × 10 −9 kg s −1 . In contrast to previous studies (e.g. in red van der Meer et al. [42], in green Daeichin et al. [51]), no dependence on the bubble size is found. We feel confident about the reliability of our shell viscosity estimates, i.e. the absence of biases potentially being introduced by the theoretical model used to solve the inverse problem, because we have proven with Fig.  3(c) that the model employed is capable of describing the (nonlinear) dynamics of single microbubbles and, most importantly, with Fig. 3(b), that it is also able to faithfully represent the overall microbubble behaviour across the whole size range with a single set of parameters, and in particular, to match the position of the experimental resonance peak, which is crucial for avoiding systematic errors. From this perspective, one can finally understand the importance of taking into account even the small systematic error R err in the sizing of the bubbles, which causes a minimal but still significant shift in the resonance peak position.
It should be noted that these values of shell viscosity do not necessarily represent the ground truth, but instead, a single-number portrayal of the complex dissipative processes occurring in the shell through the lens of the Marmottant model. Nothing prevents new rheological developments concerning the effect of dynamical parameters (e.g. the strain of the surface area) on the shell viscosity, from changing this quantitative picture.
We do not provide the result of the inverse problem regarding the dilatational modulus because it would be inconclusive. As explained in the previous section, at these acoustic driving conditions, for a wide range of bubble sizes, the elastic regime of the shell plays a marginal role and therefore the modulus can take on an arbitrary value. To correctly infer the modulus, it is necessary to probe the linear response of the bubble. Unfortunately, resolving such small oscillations exceeds the limits of optical microscopy. One must therefore resort to laser scattering techniques, which, however, compromise on the side of versatility. Notwithstanding, the range of values found in this study, as shown in Fig. 3(b), matches well the measurements carried out with the aforementioned laser scattering techniques (see e.g. Daeichin et al. [51]).

A. The shell dilatational viscosity shows large variability and no dependency on the bubble size
We believe that the variability in the bubble response is genuine and not a byproduct of the experimental methodology in light of the following arguments: (i) absence of significant arbitrary interference on the bubble dynamics from the pressure waves that reflect off from the top surface of the test chamber owing to the comparable acoustic impedances of water and the membrane material (lowdensity polyethylene, z = 1.79 × 10 6 Pa s m −1 ) and as confirmed by the pressure pulse recorded inside the test chamber with a needle hydrophone inserted through a side opening and shown in Fig. 3(a); (ii) absence of significant arbitrary interference on the bubble dynamics from pressure waves being scattered by neighbouring bubbles by reason of the extremely low dilution ratio; (iii) absence of arbitrary thermal effects affecting the shell rheological properties caused by the optical trap because of the very low localised temperature increase, estimated in the order of a few tens of millikelvin based on the work by Peterman et al. [85] and the very low electromagnetic absorption coefficient of water at the wavelength of 532 nm [82,83].
Most probably, the variability in bubble response can be rooted in the arbitrariness of the rheological properties of the shell. In support of this theory, it is noteworthy that the variability of the bubble response changes in accordance with the influence exerted by the shell properties on it, as clearly illustrated in Fig. 3(b). The rheological properties of the shell, in turn, are dependent on its microstructure. At surface tension values close to zero, a phospholipid monolayer self-assembles in a twodimensional polydomain liquid crystal with different tailgroup orientations. High variation in the size and distribution of such phospholipid domains on the bubble interface has been observed by Borden et al. [86] by using fluorescence and electron microscopy and by Kooiman et al. [87] with confocal fluorescence microscopy. The edge of the domains, where phospholipids with different tailgroups interact with each other, can be considered as defects of the bubble shell. A larger amount of defects results in more frequent interactions between phospholipids with different tailgroup orientations which, in turn, could yield a higher shell viscosity as experimentally proven by Hermans & Vermant [27] by testing the interfacial shear rheology of DPPC at different preshear values (i.e. at different levels of crystalline microstructure refinement). Based on these findings, we can conclude that the rheological properties of the shell and consequently the acoustic behaviour of the bubble are uniquely determined by the characteristic processing history of the shell. It follows that we cannot rule out the possibility that the repeated driving pressure pulses to which the bubbles in the test chamber are subjected during the test may change the original morphology of the shell (through clustering, partitioning or even shedding off of shell material). Nevertheless, the time interval between each ultrasound pulse in our experiment (a few minutes) is designed to allow for the completion of any dissolution processes and the restoration of full lipid coverage at the bubble interface (i.e. null surface stresses). Therefore, we believe that the effect of shell aging, which inevitably occurs during the operation of microbubbles in any case, does not significantly contribute to the observed variability in shell properties that is already naturally present.
The independence of shell viscosity on the bubble size found in our measurements is consistent with the unchanging shell thickness (consisting, most probably, of only a monolayer of phospholipids) and therefore shell properties with the size of the bubble. However, our shell viscosity estimates contradict previous studies that reported a smaller scatter and an upward trend in shell viscosity as the bubble equilibrium radius increases, albeit with varying rates across studies [36-40, 42, 45-47, 50, 51]. We believe that these divergent results are an artifact caused by the measurement techniques adopted. In the next sections, we present arguments supporting our hypothesis. Most of studies which found the upward trend with the bubble equilibrium radius have estimated the shell viscosity by analysing the bubble response in the frequency domain [36-38, 42, 45, 50, 51], adopting the "bub-ble spectroscopy" technique introduced by van der Meer et al. [42] or an equivalent of it for acoustic attenuation measurements [36,37]. The bubble spectroscopy method propounds to estimate the shell viscosity from the halfenergy bandwidth of the microbubble frequency response treating the microbubble as a linear harmonic oscillator. We now show how the upward trends are caused by the overestimation of bandwidth due to experimental biases. To align ourselves with these previous works and for ease of exposition from here on, we adopt the polytropic process approximation for the gas contained in the bubble [i.e., Eq. (6)].
If one assumes that the bubble undergoes small radial oscillations, Eq. (1) can be linearised and reduced to: where x(t) = (R(t) − R 0 ) /R 0 is the normalised radial excursion, ω 0 = 2πf 0 is the natural angular frequency of the system expressed as and δ is the damping coefficient consisting of three contributions: the contribution from sound reradiation the contribution from the liquid viscosity and the contribution from the shell dilatational viscosity Assuming a harmonic acoustic driving p d (t) = p a sin(ωt) and applying the Fourier transform to Eq. (12), one obtains the frequency response H(ω) of the bubbleoscillator: where ω = 2πf is the ultrasound angular frequency. If the damping δ is small enough (δ < 0.1), this can be derived from the frequency response energy |H(ω)| 2 as where FWHM H 2 is the full width of the frequency response energy at half peak response (half-energy bandwidth) and f res = f 0 (1 − 2δ 2 ) 1/2 is the peak response frequency or resonance frequency. Finally, the shell viscosity κ s can be derived from the damping δ as Exemplary frequency response energies |H(ω)| 2 for microbubbles obeying the linearised model [Eq. (17)], with equilibrium radii R 0 = 2 − 5 µm and shell viscosities κ s = 1 × 10 −9 kg s −1 and 1 × 10 −8 kg s −1 , are given in Fig. 5(a). These two values are chosen because they represent approximately the limit values of shell viscosity found in the present study. It is good to be reminded, however, that this linear approximation only holds for markedly small acoustic driving. Based on the Marmottant shell model, as early as p a = 5 kPa, the nonlinear frequency response deviates from the linear one, as shown in Fig. 5(b) for a microbubble with radius R 0 = 3.5 µm and shell viscosity κ s = 5 × 10 −9 kg s −1 . Furthermore, this threshold in p a decreases the closer the microbubble at rest is to the buckled state (i.e., in Marmottant model terms, the closer σ 0 is to zero). This suggests that the statement made in the study of van der Meer et al. [42] regarding the absence of substantial differences in the frequency response of a microbubble driven at p a = 1 kPa and one at p a = 40 kPa is inaccurate. This oversight was caused by the use of a linearised model to simulate the frequency response of the microbubble driven at p a = 40 kPa.
In an experimental investigation, the frequency response of a microbubble can be reconstructed by scanning through a range of driving frequencies [38,42,44,45,49,50] or directly acquired using impulsive driving [51]. If the response thus obtained follows the linearised model described by Eq. (17) and depicted in Fig. 6(a) in the solid line, the physical parameters are known exactly and the response is free of experimental biases, then the bubble spectroscopy method, through the use of Eqs. (18) and (19), allows the shell viscosity to be inferred correctly, as shown for κ s = 1 × 10 −9 kg s −1 and 1 × 10 −8 kg s −1 in Fig. 6(b) in solid grey and black lines, respectively. However, if the experimental parameters used to survey the microbubble frequency response are not chosen carefully, the spectroscopy method can lead to an erroneous evaluation of the shell viscosity. Large frequency steps when scanning through frequencies, short observation times of the bubble dynamics leading to spectral leakage and high driving pressure amplitudes that would negate the linear oscillator hypothesis can cause an overestimation of the half-energy bandwidth and, as a consequence, of the shell viscosity. If we assume that these inaccuracies predominantly affect the half-energy bandwidth by a systematic deviation ∆f err , the estimated shell viscosities will be impacted by an amount κ s,err = (π/2) ρ l R 3 0 ∆f err 1 − 2δ 2 −1 . For small damping values, which is usually the case for phospholipid-coated microbubbles, this expression can be approximated as: It is clear that the cubic dependency on the equilibrium bubble radius implies that even small systematic errors on the half-energy bandwidth yield significant variations in the shell viscosity across a range of bubble sizes. The impact on the shell viscosity estimate of an experimental systematic error ∆f err = 100 kHz in the half-energy bandwidth of frequency responses, as illustrated in Fig.  6(a) in the dashed line, is given in Fig. 6(b) as dashed lines. This result can explain the unphysical increase of shell viscosity with the bubble size and the lower variability compared to the present work that many studies have observed. We would also like to clarify, in view of some confusion in previous studies, that to correctly estimate the damping δ from the frequency response full width at half maximum (FWHW) the frequency response energy |H(ω)| 2 has to be used. To estimate the damping from the frequency response magnitude |H(ω)|, showcased in Fig. 6(a) in the dotted-dashed line, the full width at 1/ √ 2-maximum must be employed instead. The consequences of estimating the damping from the halfmagnitude bandwidth instead of the half-energy bandwidth can be evaluated computing from Eq. (17) the fre- , which for small values of damping can be approximated as f HM ≈ f res 1 ± √ 3δ , as opposed to the frequency at half energy f HE ≈ f res (1 ± δ). This translates into an error in the bandwidth ∆f err = 2 √ 3 − 1 f res δ and a departure from the true value of the shell viscosity equal to: The effect of this error on the shell viscosity estimate is given in Fig. 6(b) as dotted-dashed lines. Again, the error causes an upward shift and an upward trend in the shell viscosity estimate as the equilibrium radius of the bubble increases. We also emphasise that an inappropriate choice of the bubble dynamics model or inaccuracies in the physical parameter measurements can also lead to unphysical trends in the shell viscosity.
We now examine two previous studies, in which the bubble spectroscopy method was employed to infer the shell viscosity, in order to verify whether the values found are correct or are, instead, possibly a byproduct of the inaccuracies listed in the previous section. The two studies chosen are those of van der Meer et al. [42] and Daeichin et al. [51] because they stand as exemplary models for a large class of other studies which employed the bubble spectroscopy method. The values of shell viscosity they found are shown in Fig. 4, along with the values from the present study. The work of van der Meer et al. pioneered the use of bubble spectroscopy and represents the linchpin for other closely related studies (e.g. Luan et al. [44], Helfield et al. [38] and van Rooij et al. [45]). The authors reconstructed the frequency response of single microbubbles point-by-point by scanning the frequency of the driving acoustic pulse produced by a broadband ultrasound transducer. At each input frequency, the bubble dynamics was recorded by using an ultra-high-speed videomicroscopy facility, the radial time evolution extracted from it and its Fourier transform computed to finally obtain the response at the input frequency. The work of Daeichin et al. represents the last and most developed iteration of the laser-scattering measurement technique, pioneered by Dove et al. [88] and already employed by Lum et al. [49] and Supponen et al. [50]. The authors here have developed a brilliant technique to characterise in one fell swoop the whole microbubble frequency response. The microbubble was probed with a broadband photoacoustic wave generated by a nanosecond-pulsed laser illuminating an optical absorber and the resulting oscillations were detected by light scattering of a continuous laser.
As a first step, we study the scaling laws of the possible deviations of the shell viscosity values reported in the two studies from an assumed true value in order to identify the nature of the potential errors made. In view of our results reported in §IV, we assume a value κ s,true = 5 × 10 −9 kg s −1 , constant across bubble sizes, as the true shell viscosity value. Figure 7 21) and (20), that the FWHM is measured from the magnitude instead of the energy of the frequency responses and the bandwidth itself is possibly overestimated, on average, by an amount ∆f err ≈ 190 kHz. The right part of Fig. 7 The results of this endeavour are presented in Fig.  8(a) which reveals with the dotted-dashed lines that their methodology leads to an inaccurate estimate of the shell viscosity: its value gets inflated, the scatter for a given bubble size is reduced and a strong dependency on the initial radius appears. The estimated shell viscosities compare remarkably well with the values reported by the authors (green dots), endorsing the idea that the apparent dependency of the shell viscosity on the bubble equilibrium radius is indeed only the result of an experimental artefact. The methodological sources of bias occurred in the experimental procedure are diverse: (i) E 1 : the linear model was applied to bubbles behaving nonlinearly (p a = 40 kPa); (ii) E 2 : the frequency response was sampled with too coarse a resolution (∆f sam = 50 kHz); (iii) E 3 : the bubble dynamics was observed for too short a period of time (T obs = 4.2 µs) leading to spectral leakage in the frequency spectrum of its radial time evolution; (iv) E 4 : the resonance peak area of the Fourier transform of the radial time evolution of the microbubble under investigation was inappropriately considered as the measure of its response at the specified input frequency. The aftereffect of each source of bias listed above is presented in Fig. 8(b) for a bubble with shell viscosity κ s = 1 × 10 −9 kg s −1 and two different equilibrium radii R 0 = 2 µm and R 0 = 5 µm. Note that the arithmetic sum of these values does not exactly return the results shown in Fig. 8(a) because the errors combine nonlinearly. To further confirm that the cause of the departure from the prescribed shell viscosity lies in the experimental parameters adopted, the solid thick lines in Fig. 8(a) show that reducing the methodological inaccuracies, detailed above, to negligible values (i.e. using p a ≈ 1 kPa, ∆f sam < 5 kHz, T obs > 100 µs) the spectroscopy method would be capable of recovering the correct shell viscosity. Note, however, that the small oscillations resulting from such a low driving pressure would be beyond the resolving capabilities of any existing optical microscopy apparatus.
In this regard, we acknowledge the great potential of the shell investigation technique introduced in the work of Daeichin et al.. The photoacoustic impulsive driving and the detection method based on laser scattering allow one to surmount the limitations due to the use of optical microscopy that the work of van der Meer et al. ran up against and, in principle, to obtain very accurate shell viscosity and modulus estimates from linearly oscillating bubbles. However, we point out that, in the work in question, in addition to having improperly extracted the damping from the half-amplitude instead of half-energy bandwidth, the observation time used is potentially too short and the driving pressure signal might be deviating from an exact Dirac impulse. These are all factors that can compromise the shell viscosity estimation and possibly explain the reported upward trend with the equilibrium radius. Finally, we turn our attention to the works conducted by Morgan et al. [39], Doinikov et al. [40] and Tu et al. [46,47] which estimated the shell viscosity of single microbubbles using an approach very similar to ours: by finding the theoretical radius-time curve that results in the minimum mean square error from the empirical curve, recorded with a streak camera in the first two cases or with a side laser scattering technique in the second two. These works also found the shell viscosity increasing with the bubble size. Clearly, the error analysis we carried out for bubble spectroscopy cannot be applied to these studies. However, we believe that their results can also possibly suffer from other methodological oversights causing the upward trend.
Morgan et al. [39] and Doinikov et al. [40] did not report whether the theoretical model used is able to capture the behaviour of the microbubbles across a range of sizes, as we did and shown in Fig. 3(b) and, in particular, whether it can identify the correct resonance radius for the driving frequency employed. Specifically, the microbubbles were driven at a high pressure amplitude (over 100 kPa) and thus into strong nonlinear oscillations, while the theoretical model employed is based on a linear representation of the coating physics. Therefore, the resonance peak from the theoretical model is expected to be at a larger radius than that found in the experiments. This leads the optimisation algorithm to assign a low viscosity to small-diameter bubbles and vice versa, and thus to an upward trend in the shell viscosity similar to that found in their study. Tu et al. [46,47] did employ a nonlinear model, however, they solely conducted indirect measurements of bubble evolution, by converting scattering light signals. This approach raises concerns regarding the reliability of the obtained bubble radius measurements. It is important to note that any inaccuracies in the radius measurement can introduce distortions in the profile of the bubble behaviour across bubble sizes. This, in turn, can result in a mismatch between the experimental data and the theoretical model employed for estimating the shell viscosity, significantly impacting the accuracy and reliability of the inferred shell viscosity. In this regard, Tu et al. also did not compare the theoretical model used with the experimental bubble response across bubble sizes. This omission prevents the resolution of lingering doubts regarding the accuracy of the shell viscosity measurement.
To gain a better understanding of the sensitivity of the shell viscosity estimate on the radius measurement error, let us consider synthetic data of the response of bubbles across different sizes [represented by dashed lines in Fig.  9(a)], generated using a model of bubble behaviour [illustrated by solid lines in Fig. 9(a)], but with a deliberate error R err = −3% introduced in the radius measurement. If the aforementioned model is subsequently employed to infer the shell viscosity, the error induced by the inaccurate bubble sizing will propagate to the shell viscosity estimate, leading to clear trends correlating with the bubble radius, as shown in Fig. 9(b). This arises from the discrepancy between the positions of the resonance peaks in the model and the synthetic data. It is evident that if the error on the radius has an opposite sign, the trends will exhibit a reversal in their slope. To avoid these undesirable effects in our work, we carefully verified the agreement between the model and the experimental data, and thus the absence of errors on either side, prior to estimating shell viscosity, as portrayed in Fig. 3(b).

VI. CONCLUSIONS
In this work, we conducted measurements of the dilatational viscosity of the shell of ultrasound contrast agent microbubbles. By combining ultra-high-speed microscopy imaging, optical trapping and wide-field fluorescence, we obtained recordings of the response of individual microbubbles to an ultrasound driving pulse with unprecedented accuracy. We validated an advanced theoretical model of bubble dynamics based on the compressible Rayleigh-Plesset equation, and augmented by the Marmottant model to account for the shell and the Zhou model to accurately track the pressure variation in the bubble, against the experimental data collected, and found excellent agreement, both with the individual bubble responses and the overall response curve across the entire range of bubble sizes. This model was therefore employed in an optimisation routine to solve the inverse problem of determining the shell viscosity of single bubbles from their radial time evolution. The resulting shell viscosity values show a large, one order of magnitude variability, ranging approximately from κ s = 1 × 10 −9 kg s −1 to κ s = 1 × 10 −8 kg s −1 , and an average value of κ s = 5.4 × 10 −9 kg s −1 . Most interestingly, we found no dependency of shell viscosity on the bubble equilibrium radius. These results contrast with more than two decades of previous studies, which have reported less variability and, very intriguingly, a strong upward trend as the bubble equilibrium radius increases. The latter behaviour is considered to be unphysical, nevertheless, it has thus far never been either explained or proven as such. We believe that the root cause of this issue in the majority of studies lies in the high sensitivity to errors of "bubble spectroscopy", the method employed to extract shell viscosity from the bubble frequency response bandwidth. We demonstrated with analytical arguments that an overestimation of the bandwidth leads to inflated shell viscosity values, reduced dispersion for a given bubble size, and the appearance of a strong dependence on the equilibrium radius. We pinpointed several methodological inaccuracies in previous studies that lead to an overestimation of bandwidth, including: the use of linear models to describe the nonlinear bubble response, a too short observation time of the bubble response, a too low frequency sampling rate, the use of the half-magnitude bandwidth instead of the half-energy bandwidth, to name a few. We exemplified that resolving these inaccuracies allows the same trend found in our measurements to be recovered. Furthermore, our findings indicate that the unphysical trends observed in the remaining studies, where the shell viscosity was determined by fitting radius-time curves, may be attributed to errors in bubble sizing. This issue leads to a misalignment of resonance peaks between the experimental data and the bubble dynamics model used to infer the shell viscosity. As a consequence, the optimisation procedure assigns incorrect shell viscosity values to the microbubbles, resulting in trends that falsely correlate with the bubble radius. Finally, we suggested that the variability found in the shell viscosity values reflects the arbitrariness in the size and distribution of the phospholipid domains among the bubble shells. The results presented could provide useful guidelines for characterising the shell properties of microbubbles which is critical for predicting their bioeffects in medical applications.

CONFLICTS OF INTEREST
There are no conflicts to declare.